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Claims 

[d] A method for acquisition of data where the data are col- 
lected at the nodes, which are arranged as the nodes of 
the generalized Gaussian quadratures for exponentials. 

[c2] a method for acquisition of data where the data are col- 
lected at the nodes, which are arranged as the nodes of 
the generalized Gaussian quadratures for exponentials 
on intervals in dimension one, and/or rectangles and/or 
disks in dimension two, parallepipeds, and/or balls. 

[c3] The method of claim 1 where the generalized Gaussian 
quadratures for exponentials are used as a basis for 
bandlimited non-periodic data sequences 

[c4] The method of claim 1, where the accuracy of the nodes 
of the generalized Gaussian quadratures is selected and 
the nodes corresponding to this accuracy are computed. 

[c5] The method of claim 3, where for a given bandlimit of 

the bandlimited data and a given accuracy the number of 
nodes is selected to be less than the number of data 
points. 

[c6] The method of claim 3 where prolate spheroidal wave 



functions are used as the approximating basis instead of 
generalized Gaussian quadratures. 

[c7] The method of method of claim 5, where for a specified 
approximation accuracy the nodes are computed. 

[c8] The method of claim 6, where the number of nodes is 

controlled to be less than the number of the data points. 

[c9] The method of claim 7, where the number of nodes for 
the Prolate Spheroidal Wave Functions is automatically 
computed. 

[do] The method of claim 8, where the coefficients for gener- 
alized Gaussian quadratures for exponentials used as 
approximating basis are computed from the coefficients 
of the Prolate Spheroidal Wave Functions. 

[cH] The method of claim 9, where the generalized Gaussian 
quadratures are used as interpolating basis for bandlim- 
ited data. 

[d2] The method of claim 7 where it is applied to non- 
periodic data. 

[d3] The method of claim 1, where it is applied to seismic 
data with arbitrary linear moveout. 

[d4] The method of claim 4 where it is applied to seismic data 



with arbitrary linear moveout 

[d5] The method of claim 9 where the seismic data are com- 
puted at times different than the discrete time samples 
of the given seismic data set. 

[d6] The method of claim 1 applied to spatially aliased seis- 
mic data by applying a linear moveout correction parallel 
to the dominant dip of the seismic data in a local data 
window moving around the seismic data volume. 

[d7] The method of claim 4 applied to spatially aliased seis- 
mic data by applying a linear moveout correction parallel 
to the dominant dip of the seismic data in local data 
window moving around the seismic data volume 

[c18] The method of claim 7 applied to spatially aliased seis- 
mic data, wherein selecting the accuracy and the number 
of the nodes corresponding to this accuracy, applying a 
linear moveout correction parallel to the dominant dip of 
the seismic data in local data window moving around the 
seismic data volume, effectively flattening the seismic 
data, computing the coefficients of the Prolate 
Spheroidal Wave Functions, estimating the seismic data 
for this approximating basis and undoing the linear 
moveout applied to the local data window. 
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METHOD AND APPARATUS FOR EF- 
FICIENT DATA ACQUISITION AND 

INTERPOLATION 

Abstract 

An acquisition system that has a minimal number of nodes for a 
fixed bandwidth of the measured data and a given accuracy of 
measurements is disclosed. 
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Description 

METHOD AND APPARATUS FOR 
EFFICIENT DATA ACQUISITION AND 
INTERPOLATION 

Contextual Summary of the Invention 

[0001] | n a spectrum of very diverse industries, data sequences 
are recorded and users of these data wish to approximate 
the underlying functions at points other than the data 
recording points. Some examples of data acquired/ 
recorded include but are not limited to seismic data, grav- 
ity data, magnetic data, digital and scanned pictures, 
moving pictures, etc. 

[0002] Typically the data sequences are sampled at equally 

spaced nodes like in the case of digital pictures and the 
typical methods used currently for approximating/interpo- 
lating these data sequences are based on Fourier methods 
and/or on polynomials. There are two problems with us- 
ing Fourier and/or polynomials for approximating non- 
periodic data. First, the approximation at the edges of the 



data is of poorer quality than that in the middle. Second, 
the number of nodes necessary for a quality approxima- 
tion is higher than minimally possible. 
[0003] a preferred embodiment of the present invention deals 
with both problems, by providing a high quality approxi- 
mation at the edges and at the same time using the mini- 
mal number of nodes for the approximation. The teach- 
ings of the present invention may be applied in many dif- 
ferent industries and to many different types of data, 
since they enable a computationally efficient data approx- 
imation /interpolation for general types of non-periodic 
non-equally spaced band-limited data with no additional 
assumptions. 

[0004] a preferred embodiment of the present invention applies 
the general teachings of the present invention to the area 
of seismic data acquisition and data processing. 

[0005] Almost all geophysical exploration, and in particular hy- 
drocarbon exploration, includes the use of seismic meth- 
ods. Seismic exploration generally begins with a seismic 
data acquisition in an area that has been identified as 
promising for hydrocarbon exploration. Seismic data ac- 
quisition surveys use acoustic sources (generally referred 
to as "shots") as a source of seismic waves. Those seismic 



waves propagate radially through the ground in accor- 
dance with the acoustic impedance (analogous to electri- 
cal impedance in electric circuit theory) of the geologic 
layer(s) through which the waves travel. For example, in 
Figure 1, point S represents the location of such an acous- 
tic source with respect to a vertical cross-section of 
ground showing two geologic layers or "strata" 100 and 
102. Line 104 represents the direction of travel for a point 
on the radiating seismic wavefront generated by the 
source at S. Point R lies on the interface between two geo- 
logic layers having different acoustic impedances. Accord- 
ing to basic principles of physics, when a wavefront comes 
into contact with the interface between two media of dif- 
ferent impedances (i.e., an impedance "mismatch"), at least 
a portion of the wave energy is reflected back from the in- 
terface in the general direction of the source. Generally 
speaking, a portion of the wave energy will also be trans- 
mitted across the interface, although the direction of the 
transmitted wavefront may change. This change in direc- 
tion is known as diffraction. Well-understood physical 
laws, such as Snell's law, govern the velocity and direction 
of these reflected and refracted waves. Thus, because the 
velocity and direction of incident and reflected or re- 



fracted waves can be determined mathematically, it is 
possible to identify the locations of acoustic impedance 
mismatches by measuring the amount of time it takes for 
a reflected/refracted wave to reach an observation point 
on the surface. 

[0006] For example, in Figure 1, an incident wavefront travelling 
in the direction of line 104 will reach the interface between 
layer 100 and layer 102 at point R, at which point the 
impedance mismatch between layer 100 and layer 102 
causes a reflection of the wave to travel in the direction of 
line 106 to point D (the detection point) where the arrival 
time of the reflected wave can be measured. Since the ve- 
locity and direction of travel of the incident and reflected 
wave can be determined mathematically, the depth of 
point R, as measured from the surface can be determined 
from the arrival time of the reflected wave. This process is 
known as reflection seismography, and it provides infor- 
mation about the locations, shapes, and material compo- 
sitions of various geologic features. Knowledge of these 
features may be used for locating hydrocarbons or other 
mineral resources, as well as for other uses, such as de- 
termining the geologic structure of a particular site for 
civil engineering purposes, for example. Another applica- 



tion for this general type of technology is in the medical 
field, where ultrasonic acoustic waves are used in a similar 
fashion to perform medical imaging (e.g., sonograms). 
[0007] Returning now to the problem of geophysical exploration, 
however, it is customary to employ multiple seismic de- 
tectors at different points will be used in conjunction with 
a single acoustic source, to allow seismic data to be ob- 
tained over a broad area. Different types of acoustic 
sources are used in different arrangements, depending on 
the environment in question. Typically in onshore areas, 
"Vibroseis" acoustic sources are used to transmit seismic 
waves on the subsurface, which are transmitted and re- 
flected off several subterranean layers, and eventually a 
portion of the originally transmitted energy is reflected 
back towards the earth's surface and received at detectors 
(receivers), which are spaced at predetermined spatial po- 
sitions as to minimize the acquisition cost and increase 
the seismic data acquisition survey resolution. (Vibroseis 
was formerly a registered trademark of Conoco, Inc., but 
is now recognized as a generic term in the art). For off- 
shore areas, the seismic vessel performing the seismic 
data acquisition uses airguns or waterguns which gener- 
ate a significant pressure wave that propagates in the 



subsurface. The reflected seismic energy travelling back 
from the subsurface towards the ocean bottom is either 
received at receiver streamer cables towed by the seismic 
vessel or by ocean-bottom receivers placed by oil and gas 
companies. 

[0008] Seismic survey data can also be categorized by the di- 
mensionality of the data. "Two-dimensional" seismic data 
is obtained by placing the detectors in a single line. The 
information obtained in a two-dimensional survey pro- 
vides the same type of visual perspective as Figure 1, 
where the two dimensions are linear position along the 
line of detectors (horizontal) and depth (vertical, plotted 
downward). One of ordinary skill in the art will recognize 
that the depth coordinate is interchangeable with time, 
since the arrival time of a reflected wave determines the 
depth of the reflector. Mathematically speaking, two- 
dimensional seismic data is represented as two dimen- 
sional scalar field, where the scalar value represents a 
magnitude of the seismic signal received at a particular 
surface position at a particular time. 

[0009] "Three-dimensional" seismic data is obtained by arrang- 
ing the detectors over a two-dimensional area on the sur- 
face. Usually, the detectors are arranged in some form of 



grid. A set of three-dimensional seismic data is a three- 
dimensional scalar field that represents a magnitude of 
the seismic signal received at a particular surface position 
at a particular time. During seismic data acquisition, the 
data are typically recorded in digital media, and their 
sheer volume, particularly in the case of a three- 
dimensional survey, can easily exceed several terabytes 
(lxlO 12 bytes). 

[0010] After the raw seismic data is obtained, it is then processed 
to extract useful information, typically in a graphic format. 
A variety of seismic data processing algorithms have been 
developed over the years. These algorithms take into ac- 
count the seismic source and receiver positions, estimate 
the acoustic/elastic constants of the subsurface, and fi- 
nally "migrate" the data, meaning that they identify the 
proper locations of the subsurface reflectors (i.e., the geo- 
logic features that cause the reflection of seismic waves). 

[0011] An objective of the present invention is to provide an ac- 
quisition system that has a minimal number of nodes for a 
fixed bandwidth of the measured data and a given accu- 
racy of measurements. If the measured data were peri- 
odic, then the optimal distribution of nodes is that of 
equally spaced nodes, where the minimal number of 



nodes is given by the Nyquist criterion. Specifically, it is 
sufficient to have two nodes per highest wavelength (or 
wave number) that we want to measure. 

[° 012 ] Since data periodicity rarely occurs in the real world, data 
periodicity is often enforced by applying smooth windows 
on the edges of the data, thus reducing the useable por- 
tion of the data. Alternatively, one can use polynomial- 
based interpolation, which requires sampling at rates 
higher (by at least a factor of two in practical applications) 
than the Nyquist criterion. The present invention provides 
a method and apparatus for optimal data acquisition us- 
ing the generalized Gaussian nodes for band-limited ex- 
ponentials. These generalized Gaussian nodes and corre- 
sponding weights are constructed for a given bandwidth 
and for a given accuracy of measurements. It can be 
shown that for a large number of nodes, the sampling rate 
of this method for non-periodic data approaches the 
Nyquist sampling rate for the periodic data (2 nodes per 
wavelength), and is optimal for a given bandwidth and a 
given accuracy of measurements. 

[0013] Th e generalized Gaussian nodes can be computed via the 
algorithm described in the herein-incorporated reference 
BEYLKIN, Gregory. On generalized Gaussian quadratures 



for exponentials and their applications. Applied and computa- 
tional Harmonic Analysis, 2002 vol. 12, no. 3, p. 332-373. Al- 
ternatively, the generalized Gaussian nodes can be com- 
puted via the algorithm described in the herein-in- 
corporated reference XIAO, H. Prolate spheroidal wave- 
functions, quadrature and interpolation, inverse Problems, 
2001 vol. 17, p. 805-838. 

[0014] Additional weighting can be incorporated into the nodes 
and weights using the algorithm in the aforementioned 
BEYLKIN, et al. reference. 

[° 015 ] It should be noted that the nodes of the generalized 

Gaussian quadratures do not accumulate excessively near 
the end points. It should further be noted that the rate of 
accumulation reported in the aforementioned XIAO, et al. 
reference is in error. 

[0016] | n case of malfunctioning hardware or problems with po- 
sitioning of sensors, some of the desired data points may 
be absent. A preferred embodiment of the present inven- 
tion allows one to fill in the missing information, perhaps 
using a lower band limit, provided that the number of 
missing points is not excessive. 

[0017] a preferred embodiment of the present invention com- 
prises the following steps (presented in flowchart form in 



Figure 2): 

•reading from a magnetic disk or tape or RAM or any 
other storage device, or directly from the sensors a set 
of data, which may be one, two or three-dimensional, 
or may be multidimensional (block 200); 
•selecting an accuracy with which the data is to be es- 
timated and/or interpolated; typically, this accuracy 
will correspond to the accuracy of the measurements 
(block 202); 

•selecting a bandwidth within which the data is to be 
estimated and/or interpolated; typically this is guided 
by the practical requirements of the acquisition system 
(block 204); 

•computing the nodes and the weights that are mini- 
mally required to achieve the selected accuracy and 
bandwidth or retrieving pre-computed nodes and 
weights to achieve the selected accuracy and band- 
width (block 206); 

•converting the measured data into the coefficients of 
band-limited exponentials and/or the Prolate Wave 
Spheroidal Functions and/or any other bases for 
band-limited functions. The estimation is done using 
either exact or approximate prolate spheroidal wave 



functions or interpolating basis of band-limited func- 
tions in order to control the condition numbers of ma- 
trices involved (block 208); 

•using the coefficients computed at step 5, evaluate 
the interpolated measured data at any point for any 
and/or all tasks that maybe required for said data, in- 
cluding but not limited to further processing, and/or 
graphical rendering, and/or transmission, and/or 
storage (block 210). 
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Abstract 

Wo introduce now families of Gaussian-type quadratures for weighted integrals of ex- 
ponential functions and consider their applications to integration and interpolation of 
bandlimited functions. 

We use a generalization of a representation theorem due to Carat heodory to derive 
these quadratures. For each positive measure, the quadratures are parameterized by- 
eigenvalues of the Toeplitz matrix constructed from the trigonometric moments of the 
measure, l-'or a given accuracy e. selecting an eigenvalue close to e yields an approximate 
quadrature with that accuracy, 'lb compute its weights and nodes, we present a new fast 

These new quadratures can be used to approximate and integrate other esseni ially 
bandlimited functions, such as Bessel functions and prolate spheroidal wave functions. 
We also develop, for a given precision, an interpolating basis for bandlimited functions 
on an interval. 
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On Generalized Gaussian Quadratures for Exponentials and 
their Applications 

G. Beylkin 1 and L. Monzon 2 

Department of Applied Mathematics 
University of Colorado at Boulder 
526 UCB, Boulder, CO 80309-0526 

I Introduction 

In this paper we relate the Caratheudory representation of finite sequences in terms of 
exponential sums with the computation of generalized Gaussian quadratures for expo- 
nentials. Generalized Gaussian quadratures were investigated by Markov [16, 17], Krein 
[12], Karlin [11] and, more recently, by Yarvin and Rokhlin [29]. In [29] the authors intro- 
duce practical algorithms for computing the nodes and weights of generalized Gaussian 
quadratures. The resulting approximations have a number of important applications in 
a variety of fast algorithms [30], [3]. 

The Caratheodory representation theorem asserts existence and uniqueness of the 
representation of a finite sequence of complex numbers c = (ci, c 2 , . . . , c N ), c ^ 0, in the 

*° rm M 

c* = $>e^*, (1.1) 

for k = 1, 2, . . . , N and M < N, where -1 < 8j < 1 and pj > 0. Caratheodory represen- 
tation (1.1) has been the foundation for a number of algorithms for spectral estimation, 
in particular [19] is known in electrical engineering literature as the Pisarenko method. 
In this paper we develop a fast algorithm for finding M, the phases 0 = (9\ , • • • , 9m) with 
\9j\ < 1, and the weights p = (pi, • ■ • ,Pm)- Our algorithm differs from that described 
in [19], although the basic approach is similar. We achieve finite but arbitrary accuracy 
and our algorithm requires 0(N(\ogN) 2 ) operations (0(N 2 ) in a simpler version). 

One can view finding M, the phases 0, and the weights p as a nonlinear inverse 
problem for the unequally spaced discrete Fourier transform [4], [1]. It is interesting to 
note that the associated linear problem, namely the problem where M and \9,\ < 1 in 
(1.1) are given, can be arbitrarily ill conditioned. In other words, the condition number 

Research supported in paxt by DARPA grants F49620-98-1-0491 and F30602-98-1-0154 and Univer- 
sity of Virginia subcontract MDA972-00- 1-0016. 

2 Research supported in part by DARPA grant F30G02-98- 1-0154. 
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oftheVandenrioiRlt! matrix {y *'*} w= ,, M can be arbitrarily large. On the other hand, 

the nonlinear problem is well posed and we will show the I 2 — nonri estimate 

l|p|| 2 <V^||c|| 2 . (1.2) 
The main goal of this paper is to extend Caratheodory representation and use 
it to compute quadratures for integrals involving exponentials, as well as the Bessel 
and the prolate spheroidal wave functions (PSWF). These bandlimited or essentially 
bandlimited functions play a central role in many problems of signal processing and 
numerical analysis. We also consider the associated interpolation problem involving these 
functions. Specifically, we develop methodology to represent bandlimited functions on an 
interval using exponentials {e wxt '}^i L , where the bandlimit c is a positive real constant 
and t t , \t t \ < 1, are phases computed for a given c and accuracy e. With this approach, 
we are no longer limited to representing periodic functions as is the case with Fourier 

In order to consider bandlimited functions on an interval, the PSWF were intro- 
duced in [23] and [14]. Recently, a method for computing generalized Gaussian quadra- 
tures for PSWF and, as a consequence, for bandlimited exponentials, was introduced in 
[28]. In [28] the authors first construct quadratures for the I 'SWF using the fact that the 
first n of these functions form a Chebyshev system, for any n. The approach for com- 
puting generalized Caussinn quadrat ures in [2!)] relies on a variant of Newton's method 
in conjunction with a continuation procedure. As a method it is quite general but is 
computationally intensive, although in many applical ions speed is not a limitation. 

Our approach (fillers from that in [28]. We first develop a method for construct- 
ing optimal nodes and weights for integrals involving exponentials and then show that 
the same nodes and weights also provide quadratures for other essentially bandlimited 
functions, e. g. the PSWF. 

The reader familiar with Gaussian quadratures should be warned that our method- 
ology for generating quadratures is substantially different from the existing methods. The 
resulting quadrature formulas do not coincide with any other existing quadrature but, 
numerically, they are close to those in [28] . 

As a method for constructing generalized Gaussian quadratures, our results are 
limited to integrals (with a fairly arbitrary measure) involving exponentials. Our algo- 
rithm involves fini ling i igi nvalues and eigenvectors of a Toeplitz matrix constructed from 
trigonometric moments of the measure and then computing the roots on the unit cir- 
cle for appropriate eigenpolynomials. In particular, each eigenpolynomial with distinct 
roots gives rise to an identity which, for small eigenvalues, provides us with a Gaussian- 
type quadrature and also with a representation of positive definite Hermitian Toeplitz 
matrices. In these identities the size of the eigenvalue determines the accuracy of the 
quadrature formula. 

It turns out that in the case of the weight leading to PSWF. the nodes of the cor- 
responding Gaussian quadratures are zeros (appropriately scaled to the interval [—1,1]) 
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of discrete PSWF corresponding to small eigenvalues. 

As an application, we use the new quadratures to obtain efficient approximation* 
of nonperiodic bandlimited functions in terms of linear combinations of exponentials. In 
fact, we consider integrals of the form 

u(x) = jy°*> dti{t), (1.3) 

where d/j,(t) is a measure, typically d/j,(t) = w(t)dt, with w a weight function (w is a real, 
non-negative, intcgrable function with w{j)dr > 0). For a given bandlimit c> 0 and 
accuracy < > 0, our goal is to approximate u(x) on the interval '-), f using the sum 

«(*) = f>e ic ^, (1.4) 

where w k > 0 and M = M(c, e), so that 

|ti(i)-fi(i)|<e for ie[-l,l]. (1.5) 

Since it is appropriate to view (1.4) as a quadrature, we will refer to 9 k as nodes and w k 
as weights. 

In order to find u as in (1.4)-(l.o), we sample the function u in (1.3) at equally 
spaced points so that we satisfy or exceed the sampling rate dictated by the Nyquist 
criterion. The equally spaced samples of u can be viewed as the trigonometric moments 
of the (rescaled) measure in (1.3). We extend Caratheodory representation to find M, the 
nodes {9 k }^ , and the weights {vJ k }^ . We then use the same M. nodes, and weights to 
define u in (1.4) for all x e (-1, 1) and show that (1.5) holds if « in (1.3) was sufficiently 
oversampled. 

For the interpolation problem for bandlimited functions, we consider the linear 
space of functions £ c = {/ G L x [-l, 1] : f(x) = T,k& a ^ ibkX {«*} e ^ > 1^*1 < c, V/c.} . 
These functions are not necessarily periodic, in [—1, 1]. The interpolation problem amounts 
to representing, with accuracy e, the functions in £ c by a fixed set of exponentials 
{e" ' },._, ■ where M is as small as possible. We show that by finding quadrature nodes 
{**}*^i f° r exponentials with bandlimit 2c and accuracy e 2 , we in fact obtain, with ac- 
curacy r. a basis for bandlimited functions with bandlimit c. The connection between 
the generalized Gaussian quadratures for exponentials and the interpolation problem 
was first described in [28]. We use similar results to construct interpolatory bases for 
arbitrary accuracy e. 

The paper is organized as follows. We present a brief description of the Pis- 
arenko method to obtain the classical Caratheodory representation and we derive the 
estimate (1.2) in Section II. In Section III we discuss generalized Gaussian quadra- 
tures for weighted integrals and prove some of their properties for weights supported 
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inside [-5,5]. In Section IV we introduce new families of Gaussian-type quadratures. 
We develop a fast algorithm in Section V to compute the nodes and weights of these 
quadratures. We solve the approximation problem (1.3)-(1.5) in Section VI and use it in 
the next two sections to obtain quadratures and interpolating bases for bandlimited func- 
tions. We also discuss various examples to illustrate these results. Finally, conclusions 
are presented in Section IX and auxiliary material in the Appendix. 
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II Caratheodory representation 

Caratheodory representation solves the trigonometric moments problem and can be 
stated as follows (see [7, Chapter 4]), 

Theorem II. 1 Given N complex numbers c = (ci, c 2 , . . . , cw), not all zero, there ex- 
ist unique. M < A\ positive numbers p = (p\, p 2 , ■ ■ ■ ■ pit) and distinct real numbers 
0i, 0 2 , . . . , 9 M , -1< 0,- < 1, such that 

M 

c* = X>e ; ^* for * = 1,2 JV. (2.1) 

Although the theorem applies 1o all finite sequences c of complex numbers, it is 
useful in practical applications if there is a reason to seek representations of the form 
(2.1) with positive weights pj. For example, if the sequence c are the values of a covari- 
ance function, then this theorem provides the foundation for several spectral estimation 
algorithms in signal processing, e.g. the so-called Pisarenko method (see [19] for more 
details). In this paper, we are interested in the case where the sequence c contains the 
trigonometric moments of a positive weight. 

Given c, finding M, the phases 6 = {6 U ■■ ■ , 8 M ), \0j\ < 1 and the positive weights 
pcan l>e viewed as a nonlinear inverse problem for I lie unequally spared discrete Fourier 
transform [4], [1]. 

As discussed in the introduction, the problem of finding p where c, M, and \9j\ < 1 
are given, can be arbitrarily ill conditioned, hi contrast, the phases 9j in Caratheodory 
representation are related to the vector c and we have a stability estimate: 

Theorem II.2 Vectors c and p as in Theorem III satisfy 

||p|| a <V5||c|| a . 

For the proof see Appendix A. 

Grenander and Szego's proof of Caratheodory representation [7] provides a method 
to obtain M, the phases G, and the weights p. It is also the foundation for Pisarenko's 
method for spectral estimation [19]. We now outline its main points. 

II. 1 Algorithm I: Method to obtain M, 8, and p. 

1) Given c = (ci, c 2 , . . . , cat), we extend the definition of cj, to negative k as c_ k = cj 
and we define c 0 so that the (N + 1) x (N + 1) Toeplitz matrix T K of elements 
(T N ) kj = A;, lias nonnegalive eigenvalues and a.l leasl one eigenvalue is equal to 
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2) Define M as the rank of T N . By construction, we have M < N. We also say that M 

is the rank of the representation (2.1). 

3) Let T M be the top left principal submatrix of order M + 1 of TV. That is, the matrix 

T\ f has elements (c ; / ; ) (KJ M . Find the eigenvector g corresponding to the zero 
eigenvalue of T M . 

4) Construct the polynomial (eigenpolynomial) whose coefficients are the entries of the 

eigenvector q. As shown in [7, p. 58], the M roots of this eigenpolynomial are 
distinct and have absolute value one. The phases of these roots are the numbers 

5) Find the weights p by solving the Vandermonde system (2.1) for k = 1, • • ■ , M. They 

will, in addition, satisfy = i:„. 

Remark II. 1 With the extension of the sequence c k , (2.1) is valid for \k\ < N. If q = 
(f/o, • ■ • , q M ) is the eigenvector obtained in Tan 3) of Algoril Inn 1 1.1 . I hen 

M 

Y^c k+3 q k = 0, (2.2) 

k=o 

for all s, -N < s < 0. In other words, we have found an order M recurrence relation for 
the original sequence {c^}^, . 

Remark II. 2 In practice, we are interested in using Caratheodory representation if M is 
small compared with N, or more generally, if most weights are smaller than the accuracy 
sought. However, in such cases, TV has a large (numerical) null subspace that causes 
severe numerical problems in determining cq. the rank M, and the eigenvector q. 

Nevertheless, if the sequence c are the trigonometric moments of an appropriate 
weight, we will be able to modify the previous method in order to obtain the phases 9j 
in an efficient manner. In this setting, the phases and weights in Caratheodory repre- 
sentation can be thought of as the nodes and weights of a Gaussian-type quadrature for 
weighted integrals. Once the phases are obtained, Theorem II.2 assures that the compu- 
tation of the weights is a well-posed problem. In Section V.2 we present a fast algorithm 
to obtain the weights by evaluating certain polynomials at the nodes e 1 **'. 
Remark II. 3 Given any Hermitian Toeplitz matrix T, let us consider its smallest eigen- 
value . It is easy to see that Caratheodory representation implies the following rep- 
resentation of T as a sum of rank one Hermitian Toeplitz matrices, 

( T _ A (% = £ ft e to *'< , -* ) , (2.3) 
i=i 

where pj are positive and e i9 ' are the zeros of the eigenpolynomial corresponding to the 
eigenvalue A (JV) . 
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III Generalized Gaussian Quadratures for Exponen- 
tials 

IH.l Preliminaries: Chebyshev Systems 

In this section we collect some definitions and results related to Chebyshev systems. We 
follow mostly Karlin and Studden [11] (see also [12]). Readers familiar with this topic 
may skip this section. 

A family of n + 1 real-valued functions </„, .... u n defined on an interval / = [a. b 
is a Chebyshev system (T-system) if any nontrivial linear combination 

»(*) = £«*«*(*). ^ 

has at most n zeros on the interval /. This property of a T-system can be viewed as a 
generalization of the same property for polynomials. Indeed, the family {1, t,t 2 ,..., t n } 
provides the simplest example of a Chebyshev system. 

Alternatively, a T-system over [a, b] may be defined by the condition that the n+1 
order determinant is non-vanishing, 

«o(*o) «o(*i) ■■• 

tuft,) ••• mitn) 

u n (t Q ) «„(ii) ••• u n (t n ) 

whenever a < t 0 < h < • • • < t n < b. Without loss of generality, the determinant can be 
assumed positive. 

Let uo, ...,«„ be a T-system on the interval /. The moment space M n+ \ with 
respect to u 0 , . . . , u„ is defined as the set 

it)dii(t),j = 0,...,n}, (3-3) 



^ n+1 = {c=(c 0 ,...,c n )eR n+1 |c i = ^« ; 



where the measure n(t) ranges over the family of nondecreasing right-continuous functions 
of bounded variation on the interval L It can be shown that the moment space is a closed 
convex cone. We will denote the interior of the moment space Int(.M„ + i). 
Let us consider a representation of a point c = (c 0 , . . . ,c n ) e M„ + i 

cj = ]£w»j(**). j = 0,...,n, (3.4) 

where pk > 0, a < tk < 6, k = 1, . . . , m. The index 2(c) of a point c e M n +i is defined 
as (.he minimum number of points lhat arc used in the representation (3.1). where the 
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boundary points t k = a and t k = b are counted as 1/2 and the points in the interior of 
the interval a <t k < b are counted as 1. 

The representation (3.4) induces a generalized Gaussian quadrature for the in- 
tegral with the measure that defines the point c in the moment space while the index 
describes the number of nodes necessary for the quadrature. The following theorems 
(see [11] and [12]) generalize to any T-system the usual Gaussian quadratures for the 
polynomials {1, t, ■ ■ ■ , t n } . 

Theorem III.l A point c e M n +i, c^O, is a boundary point of M n+ i if and only if 
Z(c)<(n + l)/2. 

Theorem III. 2 Let u 0 , . . .,u 2m be a T-system on [a,b] and c be a boundary point of 
M.2m+i ■ Then there exist a unique representations with the index less than m + 1 which 
involves no more than m + 1 nodes. 

Theorem III. 3 Let u 0 , . . . , u 2m be a T-system on [a, 6] and c be an interior point of 
Mim+L- Then there exist at least two represe.ii.lulioies iriJh I lie re.de.x m + 1/2 (with m + 1 
terms). Both of them have m + 1 nodes one of which is the end pond of the interval. 

If the functions u Q ,.. .,u n are periodic on the interval / and satisfy (3.2), then 
they define a periodic T-system. A periodic T-system always involves an odd number 
of functions. Indeed, since the system is defined on a circle, the equally spaced values 
(t 0 , •••,t n ) can be continuously rotated into (/,„, • • • , /. 0 ). If the number of functions were 
even, such rotation would change the sign of the determinant in (3.2) and, due to the 
continuous dependence on (to, • • • , t n ). would force the determinant to vanish at some 
intermediate point. 

For periodic systems holds ([11], [12]): 

Theorem III. 4 Let Uo, ■ ■ ■ , «a m be a periodic T-system on [—1, 1] and c be an interior 
point of Mim+i. Then for each point t 0 , -1 <t 0 <1, there exist a unique representation 
with the index 1(c) = m + 1 (with m + 1 terms) involving to as a node. 

III. 2 Generalized Gaussian Quadratures for Exponentials 

Let us consider a family of real periodic functions 

l,cos(7rt),sin(7rt), . . .,cos(7rm*),sin(7rmt) (3.5) 

on the interval [—1,1]. We treat the boundary points -1 and 1 as identical so that (3.5) 
is defined on the circle. The system of functions in (3.5) is a periodic Chebyshev system 
(T-system) in [-1, 1] (see [11], [12] and Section III.l for a brief summary) . 

We also consider this Family on a proper snliinlerval of [-1, 1]. On any subinterval 
[a, b] C [-1, 1] the family in (3.5) is a T-system. 
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Let us consider the moments of the measure ^(r)dr, where w(t) is a weight 
function, 

ak = j w(r)cos(7rikT)dr, Jb = 0,l,...,m, (3.6) 
and i 

A = y w(T)rin(jr*r)dr, fc = l,...,m. (3.7) 

We also consider complex-valued moments, 

lk = a h + \3 k = y ^ wMef^dT, A; = 1, ... ,m. (3.8) 

Let A4 2 m+t be the moment space and Int(A4 2m+ i) its interior, as defined in Sec- 
tion III.l. We have 

Theorem III.5 ([11], VI, Sec.4) For the periodic T-system (3.5) a point 
c = {ay,, a u Bi, . . . , a rn , 3 m } is a point of the moment space Mim+i if and only if the 
Toeplitz matrix {^k'-k}k,k'=u,..,m is non-negative definite. 

Furthermore, c e Int(«M2m+i) if end only if the Toeplitz matrix {jk-k'}k,k'=o,...,m 
is positive definite. 

We also have 

Theorem III.6 ([11], VI, Sec.2) A point c - {c* 0 , «i> Pi, ••• ,«m,An} is a boundary 
point of the moment space M' 2m +i if and only if there is a unique representation 

Ti = |>e™'*, (3.9) 

where m' <m and -1 < 8j < 1. 

If c e Int(A^2m+i)> then for each tq G [—1, 1] there exists a unique representation 
with m + 1 nodes including r 0 as a node, that is 

7* = f>e^ W"" 0 *, (3.10) 

where -1 < 9j < 1. 

Lei us consider weigh! s supported in a subinterval of [-1, 1]. We then prove 
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Theorem III. 7 Let to be a wt i > < ipp rtt I in some interval I = [a,b], I C [—1,1]. 
Then, there exists a uniqai representation 

j'wit^dt^w^ + ^i-lf for \k\<m, (3.11) 

where a < 9j < b and uij > 0 for j = 0, . . . , m. 

Moreover, if I = [—0,0], where 0 < a < 1/2, then 



(3.12) 



°- 2 + (2 + VS) ra + (2-V3)"*' 

Proof 

Let us start by considering the periodic T-system (3.5) and 
c = {o-o, a-\, 0-\, . . . , a m , f3 m } the point in the moment space obtained from 
d/j,(t) = ui(t)dt. It is easy to show that the Toeplitz matrix {"fj-k}, obtained for the 
moments in (3.8), is positive definite (see (4.7)). Thus, Theorem III.5 and (3.10) with 
T 0 = 1 imply a unique representation (3.11) with —1 < 0, < 1 and k = 0, •••,»». Since 
the weights are real, (3.11) also holds for k = -m, We want to show that, in 
fact, a<9j< b. 

Since the weight lo is supported in [a, b] C [—1, 1], we can also consider [a, 1] as its 
interval of definition. We note that the functions in (3.5) form a Chebyshev system on 
this interval (in fact, on any subinterval of [—1, 1]). Using Theorem III.3 we construct 
a representation which includes the boundary point 1 of the interval [a. 1] as a node. 
However, this representation also holds on [—1,1] where (3.11) guarantees uniqueness. 
Thus, to avoid a contradiction, we conclude that a < 9j < 1 in (3.11). Similarly, let us 
consider the weight u in [-1, 6]. By the same argument we obtain that -1 < 9j < b (the 
points -1 and 1 are identical on [-1,1]). Therefore, we conclude that a < 9j < b and 
(3.11) is established. 

Let us now consider a periodic trigonometric polynomial 

v m (t) = T m (l - cos(rf)), (3.13) 

where T m is the Chebyshev polynomial of degree m. Since the degree of v m (t) does not 
exceed m, we have from (3.11) 

J ' v m (t) u(t) dt = f>U%) + (3-14) 



where -1/2 < 9j < 1/2. 
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Lei uy compute v m (l) = 7„,(2). Using the three-term recursion for the Chebyshev 
polynomials, we obtain 

v m {l) = ((2 + V3) m + (2 - V3) m )/2. (3.15) 
Since in the interval [-1/2, 1/2] the absolute value of v m does not exceed one, 

0*1^(1) = v m (t)u,(t)dt-'pu j v m (O j )^ < j\ u{t)dt + , pu J . (3.16) 

Setting k = 0 in (3.11), we obtain 

f> = J\(t)dt- (3.17) 

and, combining with (3.16), we arrive at (3.12). □ 

Remark III.l Since the numerator in (3.12) remains bounded and the denominator 
grows exponentially last with in. the coellicient u.' 0 is very small even for m of moderate 
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IV A new family of Gaussian-type quadratures 

Let us consider the trigonometric moments of a weight w(r), 

tk = f Z mTk w{T)dT. (4.1) 

In our approach, it is essential to consider only weights supported inside -!, : ,'\ Only 
then can the moments t. k be viewed as values of a properly sampled baiidlimited function 
(see (6.2) and (6.3)). 

In this section we start by using Caratheodory representation and Theorem III. 7 
from the previous seel. ion, l.o eonsl mel I. wo difl'erenl Gaussian quadratures for integrals 
with weight w. These quadratures are exact for trigonometric polynomials of appropriate 
degree. 

We then generalize these types of quadratures furl her and develop a new family 
of Gaussian-type quadratures. This family of quadrature formulas is parameterized by 
the eigenvalues of the Toeplilz matrix 

T = {*,_*}„<*,!<"■ (4-2) 
Among these now quadrat. ure formulas, only I hose corresponding lo eigenvalues of 
small size are of practical interest. In fact, the size of the eigenvalue determines the error 
of the quadrature formula. To compute the weights and nodes of these quadratures, we 
develop a new algorithm which may be viewed as a (major) modification of Algorithm 
II. 1. The new algorithm is described in Section V. The main results of this section are 
gathered in Theorem IV. 1. 

We start by using Theorem III. 7 to write 

t» = ^« i e h **+w 0 (-l)* for \k\<N, (4.3) 

for unique positive weights uij and phases 4>j in (—1, !)• Then, for any A(z) = Y^\k\<N a k zk 
in An, the space of Laurent polynomials of degree at most N, we have 

f A(e^)w(r)dr = £ a k t k = jr M] A(e^) + uj 0 A(-1), (4.4) 

J - L \k\<N j=l 

for unique positive weights uij and nodes e"^' . 

Alternatively, using Caratheodory representation (2.1) applied to the sequence 
c k = t k ,l < k < N, 

f_A{^)w{T)dr = ^A^ + ito-c^y^A^ndT 
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= ^ Pj A{^) + \W\j\e™)dT, (4.5) 

where Cq = YljLiPj an d W* e '} are the roots of the eigenpolynomial corresponding to 
the smallest eigenvalue X^"> of T. 

Note thai (4.5) is again valid for all A(z) in A. v and that the positive weights />, 
and phases Oj in (—1, 1] are unique. 

Thus, wo have two different quadratures! that, may not coincide. However, by 
considering w(t) supported inside (-}„\), (3.12) implies that «: 0 in (4.4) decreases ex- 
ponentially last with ,V and. since min w{t) = 0 for r < 1, we have 

lim A^ = 0 (4.6) 

as shown in [7, page 65]. In consequence, for large N, the difference between these two 
quadratures can be made smaller than the accuracy sought. 

A similar reasoning could be applied to other small eigenvalues of T provided we 
can generalize (4.5) to other eigenvalues and roots of the corresponding eigenpolynomials. 
For that purpose, we first describe some properties of these eigenpolynomials. 

IV. 1 Toeplitz matrices for trigonometric moments 

We summarize in this section properties of eigenpolynomials of the Toeplitz matrix T 
with entries i =0 ... N . The matrix T is self-adjoint and positive definite since, for 

(Tx,x)= J2 fjP x {^ T )\ 2 '<r)dr (4.7) 

where (sc, y) = J2 k x kVk is the usual inner product of two vectors and P x {z) = ^2 k x k z k . 

More generally, for all x, y £ C N+1 , T induces a weighted inner product for 
trigonometric polynomials, 

{Tx, y) = £ P x (e^)P^)w(r)dT. (4.8) 

Since T is positive definite, there exist an orthonormal basis {w ( *'}fc=u,...,jv of eigenvectors 
of T corresponding to eigenvalues A (0) > A (,) > ■ - > A (,v) > 0. The corresponding 
eigenpolynomials V {k) [z) = X), vf'z' satisfy 

J vW(e i * t )VW(e^)dt = 2$ kl (4.9) 



Gnussinn qumimuircs for oxponciitiak 



14 



and, because of (4.8), 

J' V^(^)VW{^)w{t)dt = 5 kl \W. (4.10) 
For a vector x = (x 0 , ■ • • , x N ), let us define the reciprocal vector of a: as 
x„ = (xW, --^xo) 

and. similarly, for a trigonometric polynomial P(z). the reciprocal polynomial of P as 

P i (z) = P(z-i) = P(z^j. (4.11) 

Since T is Toeplitz Hermitian, we have 

Txi, = Ax*, if Tx = \x. 

In particular, if A is a simple eigenvalue, its corresponding eigenspace can be gener- 
ated by a self-reciprocal eigenvector x, i. c, x, = x and the associated (self-reciprocal) 
eigenpolynomial will have roots in pairs {7, 7 -1 }. 

IV. 2 Gaussian-type quadratures on the unit circle 

In this section we present the main results of the paper. We derive new Gaussian- 
type quadratures valid for any eigenvalue of the matrix T rather than just the smallest 
eigenvalue A iV . These quadratures allow us to select the desired accuracy and thus, to 
construct accuracy dependent families of quadratures. 

The nodes of the quadrature in (4.5) are the roots of the eigenpolynomial corre- 
sponding to the least eigenvalue of T and, because of Caratheodory representation, we 
know that these roots are on the unit circle and that the weights are positive numbers. In 
our generalization, t his si andard property for the nodes and weights is no longer inforced. 
However, we will show that for nodes on the unit circle, the corresponding weights are 
real. Moreover, in all examples we have examined, for all small eigenvalues A of T, their 
negative weights are associated with the nodes outside the support of the weight and are 
comparable in size with A. We believe this property to hold for a wide variety of weights. 

We prove the following 

Theorem IV. 1 Assume that the eigenpolynomial T ''•"•'(•:) com sponding to the eigen- 
value A (s) ofT has distinct, nonzero roots {7j}£Li- Then, there exist numbers {wj}j =l 
such that 

i) For all Laurent polynomials P(z) of degree at most N, 

jf ' P(^)w(t) dt = "£ WiPin,) + A«± P(e^) dt. (4.12) 
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ii) For each root 7^ with fo] = 1, co/ - /« «/*'«</ ' eight w k is a real number and 

m = £ \L k (e™)fw(t) dt - AWi £ |^(e irt )| 2 dt, (4.13) 

tuyere 

^•'- (V^'g-Tj '" 4) 
is Uu: liugnmiji: polynomial aHHwitdwl with l,fi<: fool ", k . 

iii) //A is a s*m/>/t x ' Ihtn for k = I , • • • , A ; /At weight w k is nonztro and 

Jt4 AO -AW ' t4 - i5j 

wfcere vj 0 (z) = ^(z" 1 ) is the reciprocal polynomial ofV®{z). 
In particular, for each 7^ with \j k \ = 1, 

L = T (416) 
ok^'-Aw (4 - 16) 

fie 

iv) // A (s) is a simple eigenvalue and all roots y k are on the unit circle, then the set 

{ w k}k=i contains exactly s positive numbers and N — s negative numbers. 

In particular, if s = 0 or s = N then all w k are negative or positive, respectively. 

Remark IV. 1 Our approach to obtain Gaussian quadratures does not use Szego poly- 
nomials and is therefore substantially different than the one in [10]. We briefly explain 
the approach in [10]. Note that (4.9) and (4.10) show that the polynomials {V^(z)} are 
orthogonal with respect to both the usual inner product for trigonometric polynomials 
and the weighted inner product with weight w(t). We can also construct Szego polynomi- 
als {pfc(z)} orthogonal with respect to w(t) and such that each p k (z) has precise degree 
k [25]. For any k, the roots of p k (z) are all in \z\ < 1 [7]. 

Szego polynomials and I heir reciprocals induce para.-orl hogona.l polynomials [10], 

B n (z)=p n (z)+t n z n ( Pn )*(z), 

where £„ are complex constants, |£„ = 1. The roots of B n (z) are on the unit circle and 
can he used as (.ho nodes for Gaussian quadrat urns with respect to the weight «:(*). 

Under appropriate assumptions to guarantee uniqueness, the quadratures in II) 
should coincide with those obtained in Theorem 1 1 1.7. 
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In contrast to these exact quadratures, in Theorem IV.l we derive a new family 
of Gaussian-type quadratures where for each eigenvalue of the Toeplitz matrix (4.2) 
there is a corresponding quadrature formula. Even for the smallest eigenvalue A^'', the 
quadratures in (4.12) and in III. 7 are different because of the extra integral term in (4.12) . 
The size of this extra term is controlled by the size of the corresponding eigenvalue and, 
thus, it is never exactly zero. However, in applications, this extra term can be made as 
small as desired via oversampling (see (6.1)-(6.3) and note that (4.6) is valid for all small 
eigenvalues, not just the smallest [7, page 65]). 

Remark IV. 2 For an eigenvalue A (s) of small size, the integral term on the right hand 
side of (4.12) can be neglected. This is the case of practical interest. 

Remark IV. 3 Even though the weights w k could be complex valued, an important 
consequence of Theorem IV.l is that in many important cases w k are, in fact, real. 

Remark IV.4 We have observed (see Table 1) that the weights w k corresponding to 
nodes outside the support of the weight w(i) are small, negative, and roughly of the size 
of the eigenvalue A^. Although we now present a heuristic explanation of this behavior, 
we do not know if a proof can be obtained along this path. 
Let us split the sum in (4.16), 

w k f , AW-AW ^ AW -AW [ '> 

If A w is in the range of the exponential decay of the eigenvalues, the first term in 
(4.17) turns out to bo much smaller than tiro second term which is approximately 

A^y E I^W- 

(:AM>A«> 

For 74 outside the support of the measure we have observed (Figures 2, 3, and 
5-8) that 

E i^ (i) (7*)i 2 

I:A(')>AW 

is a constant of moderate size. 

Thus, the second term in (4.17) is 0(1/A (s) ) and the weight is indeed negative 
and roughly of the size of the eigenvalue. 

Remark IV. 5 For the weight with value 1 in {-\, i) and 0 otherwise, the eigenpoly- 
nomials are the Discrete PSWF. For these functions we know that all eigenvalues are 
simple and that, all eigenpolynomial roots are on tlie unit circle >2\ 
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Corollary IV. 1 Under the assumptions of Theorem IV. 1, it follows that the Toeplitz 
matrix T in (4.2) has the follower) , , , , , la lion as a sum of rank one Toeplitz matrices, 



wlwn: AW, Wj, andjj are as in (4.12). 

This corollary should be compared with Remark II. 3 noting that, in the corollary, 
A (s) is not necessarily the least eigenvalue of T. 
Proof of Theorem IV. 1 

1. For x = {x 0 , • • ■ , x N ) e C N+1 let us define 

Ytt-TXi+LZ 1 if N = 2L 



A x (z) = 



El-L+^i+L-iZ 1 if N = 2L-l. 



The values of A x on the unit circle have a phase shirt respect of those for l' x . In fact , 
depending on the parity of N, A x {er') is either P x (e int )e- ML or P a! (e i,ri ) e ~ i,rt(i ~ 1) - 
Hence, (4.8) holds replacing P x by A x and then (4.9) (4.10) also hold for the shifted 
eigenpolynomials. 

We prove the Theorem for the case N = 2L. (The case N = 2L — 1 is similar.) For 
this case, using the same notation V (fc) for the shifted eigenpolynomials. we have 

i=-L 

2. Since {7^} are distinct, we define {ujj}f =l as the unique solution of the Vandermonde 

system 

J2 -yJ k Wj = j ' e-' mtk w(t)dt for k = 1, ■ • • , N. (4.18) 

3. Let P e A N , then z N P(z) it is a polynomial of at most degree 2JV, and since 2 l V (s >(jz) 

it is a polynomial of degree N, by Euclidean division, there exist polynomials q(z) 
and r(z) of degrees at most N and N — 1 such that. 

z N P(z) = W\z) q {z)+r{z). 
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P(z)=V^(z)Q(z) + R(z), (4.19) 
where Q(z) £ A L and R(z) has the form R(z) = J2 k =i r k z ~ k and hence 
j R(e i7rt )dt = Q. 

Using the fact that {V^}fL 0 is a basis of A L , we write 
Q(^) = f>V'«(e«), 
where d ( are some complex coefficients. 

Using (4.10) and (4.18), we multiply both sides of (4.19) by w(t) and integrate to 
obtain 

f P(e^)w{t)dt, = J2 J i f V ( ' ) (e i " t )VM(e«)w(t)dt+ f R(e mt )w(t)dt 
J-i l=0 J-i J-i 

Now, (4.19) implies that the last sum equals WjP(jj). To find the constant 
d s , we integrate both sides of (4.19) and use (4.9) to obtain 

f P(e lwt )dt = J2d, f VW(j nt )VW(e i « t )dt = 2d s , 

and thus (4.12). 

4. Let us assume that the node n / k has unit norm, |-yjb I = 1> a "d l et P( z ) = L k (z) (L s k )+(z). 

We have P(j T ) = S rk and since P € A N - U (4.12) implies 

\L k {^)\ 2 w{t)dt = X^\ y' |Li(ef rt )| a dt + «, 4 . 

Clearly w k is real. 

5. We now show that for 1 < k, j < N. 
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and thus, considering k = j, (4.15) follows. Note that we need A (s) to be simple to 
guarantee A (i) - A< s > # 0,/ # s in (4.20). 

If we view the left hand side of (4.20) as the entries A kj of a matrix A and let B 
be the matrix of entries 

B, k = V®(>y k ) where 0 < I < N, I ± s and 1 < k < N, (4.21) 
we can prove (4.20) by showing that B A = B and that B is non-singular. 
For the latter claim, we simply check that the columns of B are linearly indepen- 
dent. Indeed, let a h I s, be constants such that 

J2aiV^( 7k )=0 for k=l,...,N. 

It follows that the polynomial P(z) = ^aiV {!) (z) 6 h T , has the N = 2L distinct 

roots 7fc. Since P and V M have the same degree and the same N distinct roots, 
P(z) = c V^(z), for some constant c. By (4.9), V^(z) is orthogonal to all the 
other eigenpolynomials and so a t = 0. 

To show that BA = B,we first substitute P(z) = V«\z)v} m) (z) in (4.12) to 
obtain 

+ E»^ (i, (0;)K <m) (7i). 
Using (4.9)-(4.10), we rewrite the previous equation as 

<WA (!) - AM) = £ w^Hj^Hji) (4-22) 

and thus, 

iBA U - 

= V^{ ln ) = B mn . 
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6. To prove the last assertion of the theorem, we consider (4.20) when all 7* have unit 
norm and thus all w k are real. In this case, 

and we can rewrite (1.20) as a matrix identity 

B*YB = W, (4.23) 
where B is the invertible matrix defined in (4.21), B* is its adjoint and T and 
W are diagonal matrices with real entries j ^(i) * ^( s ) } an< ^ {^"} ' 
respectively. 

Using Sylvester's law of inertia [9, Theorem 4.5.8], (4.23) implies that V and W 
have the same inertia, thai, is, the same number of positive, negative, and zero 
eigenvalues. The result follows because we assumed A <s) to be simple and then 

\W > > AC'" 1 ' > AW > AC s+, > > > AW. □ 



Techniques similar to those used in the proof of Theorem IV. 1 allow us to derive 
several results for eigenpolynomials corresponding to multiple eigenvalues or for the case 
where their roots lie outside the unit circle. Here we limit our attention to the case of 
simple eigenvalues or eigenpolynomials with all roots on the unit circle. 

W. Trench [26] has shown that both the multiplicity of the eigenvalues and the 
number of the eigenpolynomial zeros outside of the unit circle depend on the oscillations 
of the weight function w(t). We state two of the results in [26] for T as in (4.2), 

Theorem IV. 2 ([26], Theorem 2.1) If X is an eigenvalue ofT with multiplicity m, 
then w(t) — A changes sign al least '2 m — 1 limes in (—1, 1). 

Theorem IV.3 ([26], Theorem 3.1) Letu(z) be a self-reciprocal eigenpolynomial cor- 
responding to the eigenvalue A ofT. If u(z) has 2m (m > 1) zeros that are not on the 
imil circle, then w(t) — A changes sign at least 2m + 1 times in (—1, 1). 



IV.3 Examples 

We consider three examples with diU'erem weigh is and conslrucl I he appropriate quadra- 
tures. 

Example 1. First we consider the weight 
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For this weight the eigenpolynomials V^(e mt ) of the AT + 1 x N + 1 Toeplitz matrix 
T are the discrete PSWF [22], Thus the eigeiipolyiiomial V^(e M ) has all of its zeros 
on the unit circle. Moreover, it has exactly I zeros for t in the interval (—a, a) and N 
zeros for t in [-1,1]. In this example we have selected N = 97, o = 1/6, c = 157T. 
We then construct the matrix T and compute the eigenpolynomial corresponding to the 
eigenvalue 

A (30) = 9.77306136381891632828 • 10" 16 . (4.25) 
The eigenpolynomial V' (3U ' (e 17rt ) is shown in Figures 2 and 3. Locations of the zeros on the 
unit circle are displayed in Figure 4. We then use the quadrature formula corresponding 
to this eigenvalue and tabulate the weights in Table 1. Note that the weights for nodes 
inside the interval [-1/6, 1/6] are positive, and those for nodes outside this interval are 
negative and roughly of the size of A (30) . 
Example 2. We consider the weight 

In this example we have selected JV = 61, a = 1/4, c = 157T. We then construct the 
matrix T ami compute I ho eigmipolyiiominl corresponding to the eigenvalue 

A (28) = 1.11598931688523706280 • 10" 14 . (4.27) 

The eigenpolynomial V (28) (e ,,rt ) is shown in Figures 5 and 6 . 
Example 3. We consider a non-symmetric weight 

W ^ ~ { 0 ~ elsewhere. ~ (4.28) 

In this example we have selected JV = 61, a = 1/4, c = 157T. We then construct the 
matrix T and compute the eigenpolynomial corresponding to the eigenvalue 

A< 28) = 4.68165338379692121389 ■ 10" 15 . (4.29) 

The eigenpolynomial T /' 28 '(e' 1rt ) is shown in Figures 7 and 8. Although we do not have a 
proof at the moment, it appears that there is a class of weights for which eigenpolynomials 
corresponding to small eigenvalues mimic the behavior of the discrete PSWF with respect 
to locations of zeros. In Example 3 we know that all zeros are on the unit circle due to 
Theorems IV.2 and IV.3. 

In Table 2 we illustrate the performance of quadratures for different bandlimits 
c. This table should be compared with Table 1 in [28]. The performance of both sets 
of quadratures is very similar. Yet these quadratures are quite different as can be seen 
by comparing Table 3 with Table 5 in [28]. Although the accuracy is almost identical, 
approximately 10~ 7 , the positions of nodes and weights differ by approximately I0" 3 - 
-10" 4 . 
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Figure 1: Decay of the eigenvalues of the matrix T in Examples 1-3. The scale of 
the vertical axes is logarithmic (log 10 ), whereas the horizontal axes displays indices of 
eigenvalues. We note the exponential rate of decay. The flat portion of the graph for 
large indices is due to the limited precision of our computations. Thus, these graphs 
also illustrate the practical difficulty of finding the eigenpolynomial corresponding to the 
smallest eigenvalue. 



GuusHiftti qundriiUircs for oxponciitiak 



23 




Figure 2: Modified eigenpolynomial e _1 '% V^^e"") on the interval [—1,1], where 
N = 97 and l/C 30 ) (e iTt ) is the eigenpolynoiniai corresponding to the eigenvalue A (30 ' 
in Example 1. The phase factor e~ mtN / 2 i s introduced to make this function real. 




-1-10 

Figure 3: The same function as in Figure 2 on the interval [—1/6, 1/6]. 
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Figure 4: Location of the zeros on the unit circle for the eigenpolynomial V' (30) in Exam- 
ple 1. 
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Figure 6: The same function of Figure 5 on the interval [—1/4. 1/4]. 
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Figure 7: Modified eigenpolynomial (see Figure 2) on the interval [-1,1] corresponding 
to the eigenvalue A t28) in Example 3. 







'"-0.2 -0.1 

-2-10" 6 

-4-10" 5 

-6-10" 6 

-8-10" 5 





Figure 8: The same function of Figure 7 on the interval [-1/4, 1/4]. 
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# 


weights 


1 


-1.0328 • 10" 17 


2 


-1.0328 • 10" 17 


3 


-1.0329 • 10- 17 


33 


-1.3518 • 10~ 17 


34 


-1.6030 • 10- 17 


35 


0.00580295532842819966 


36 


0.01310603337477264417 


37 


0.01959211245475268191 


38 


0.02506789313597245367 


39 


0.02954323947353217723 


40 


0.03313334531810570720 


41 


0.03598544514201341779 


42 


0.03823923547752508920 


43 


0.04001188663952018400 


44 


0.04139574827622469674 


45 


0.04246105337417774134 


46 


0.04325984471286061543 




0.04382960375644760677 


48 


0.04419611220330997984 


49 


0.04437540133235668283 



# 


weights 


50 


0.04437549133235668283 


51 


0.04419611220330997984 


52 


0.04382960375644760677 


53 


0.04325984471286061543 


54 


0.04246105337417774134 


55 


0.04139574827622469674 


56 


0.04001188663952018400 




f ) i ) 3 n ') 3 9 ') 3 5 4 7 7 5 2 5 0 8 9 2 0 


58 


110 359 J4514J(lli4177<» 


59 


0.03313334531810570720 




0.02954323947353217723 


61 


0.02506789313597245367 


62 


0.01959211245475268191 


63 


0.01310603337477264417 


64 


0.00580295532842819966 


65 


-1.6030 10" 17 


60 


-1.3518 • 10" 17 


90 


-1.0329 • 10" 17 


97 


-1.0328 • 10- 17 



Table 1: Table of weights for the quadrature formula with A (30) in Example 1. The weight 
#1 corresponds to the node 7, = -1 (see Figure 4). 



c 


# of nodes 


Maximum Krrors 


20 


13 


1.2 - 1()- 7 


50 


24 


1.1 • 10~ 7 


100 


41 


1.6- 10- 7 


200 


74 


1.8 • 10" 7 


500 


171 


1.4 • 10" 7 


1000 


331 


2.4- 10" 7 


2000 


651 


1.2 • 10" 7 


1000 


1288 


3.7 -10" 7 



Table 2: Quadrature performance for varying band limits. 
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Node 


Weight 


-0.99041609489889 


2.42209284787E- 


-02 


-0.95238829377394 


5.04152570050E- 


-02 


-0.89243677566550 


6.82109308489E- 


-02 


-0.81807124037876 


7.96841731718E- 


-02 


-0.73438712699465 


8.71710040243E- 


-02 


-0.64454148960251 


9.22000859355E- 


-02 


-0.55050369342444 


9.56668891250E- 


-02 


-0.45355265507507 


9.80920675810E- 


-02 


-0.35456254990620 


9.97843340729E- 


-02 


-0.25416536256280 


1.00930070892E- 


-01 


-0.15284664158549 


1.01641529848E- 


-01 


-0.05100535080412 


1.01982696564E- 


-01 


0.05100535080412 


1.01982696564E- 


-01 


0.15284664158549 


1 .0161 1 5298 18K- 


-01 


0.25416536256280 


1.00930070892E- 


-01 


0.35456254990620 


9.97843340729E- 


-02 


0.45355265507507 


9.80920075810K- 


-02 


0.55050369342444 


9.56668891250K- 


-02 


0.64454148960251 


9.22000859355E- 


-02 


0.73438712699465 


8.7 1 7100-10243I-:- 


-02 


0.81807124037876 


7.9684 17317181-:- 


-02 


0.89243677566550 


6.821 09308 1891-:- 


-02 


0.95238829377394 


5.04152570050E- 


-02 


0.99041609489889 


2.42209284787E- 


-02 



Table 3: Quadrature nodes for exponentials with maximum bandlimit c = 50. The 
maximum error is « 1.1 • 10 -7 . 
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V A new algorithm for Caratheodory representation 
V.l Algorithm 2 

We now describe an algorithm for computing quadratures via a Caratheodory-type ap- 
proach based on Theorem IV. 1. It is easy to see that, although there are similarities 
with Pisarenko's method, the algorithms are substantially different. We plan to address 
implications for signal processing in a separate paper. 

1) Civen / ( ., the trigonometric moments of a measure, we construct the {N+l) x (iV + 1) 

Toeplitz matrix T N with elements (T N ) kj = tj_ k . This matrix is positive definite 
and has a large number of small eigenvalues. 

2) For a given accuracy e, we compute the inverse of the Toeplitz matrix T N — el. For 

a selfadjoint Toeplitz matrix, it is sufficient to solve (TV — el)x 0 = e 0 , where 
e 0 = (1. 0, . . . , 0) ( . After x 0 is found, we use the Gohberg-Semeneul representation 
of the inverse of the Toeplitz matrix [6] (see also [5] for a modern perspective) in 
order to apply it to a vector. If e is too close to an eigenvalue of T N , it might be 
necessary to slighty modify the value of e and repeat this step. 
This step requires? 0(N 2 ') opera! ions if we use the Levinson algorithm. However, 
we know how to build a stable 0(N(\ogN) 2 ) algorithm for this purpose which we 
will present elsewhere. 

3) Using the power method for (TV - we find an eigenvalue A (?) close to e and 

the corresponding eigenvector q. This step requires G(N log N) operations due to 
the Gohberg-Semeneul representation of the inverse. 

4) Next, compute all zeros on the unit circle of the eigenpolynomial corresponding to 

the eigenvector q. 

This requires 0(N log N) operations since we use the unequally spaced fast Fourier 
transform ([4], [1]) to evaluate the trigonometric polynomial on the unit circle. We 
pick out the zeros within the support of the measure and denote their number by 
M. 

5) Using the algorithm described below, we find the weights p by solving the Vander- 

moroks system for all nudes (including I hose outside the support). This algorilhm 
takes 0(N log N) operations. 

V.2 Solving Vandermonde systems using polynomial evaluation 

To obtain the weights in step 5) of Algorithm V.l, we need to solve a M X M Vandermonde 
system with nodes on the unit circle. In this section we discuss an algorithm to obtain 
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the solution by evaluating certain polynomials on the Vaiidennonde matrix nodes. The 
algorithm can be derived from more general results [8, 15, 20]. Here we give a simpler 
presentation adapted to our particular application. Note that general algorithms to solve 
Vandermonde systems are unstable, unless there is a particular arrangement of the nodes. 
Let {71, . . . , 7m} be distinct complex numbers and define 



Since the nodes {7,.} are distinct, det V / 0 and thus given b = (6 0 , 
a unique p — (pi, ■ • • , p M )* such that 



, bM-iY, there is 
(5.1) 



ceo = no* 

k=l k=0 

then, for any polynomial P of degree at most M — 1, 



P{Z) ^ P(7r) 



(5.2) 



(5.3) 



Now choose P to be the unique polynomial with P(7 r ) = Pr Q'{jr) for 1 < r < M, and 
let B(z) = ^2 b * zk - Substituting in (5.3), 



= ^ p(r)~>r z k + higher powers 



= z M - L P(z- L ) and 
= z M Q(z-% 
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P(z) = Q(z)B(z)+z M Q(z)(-..) 

and thus the coefficients of P correspond to the first M coefficients of Q(z)D(z) 
As a result, we have the following algorithm: 



(5.4) 



V.3 Algorithm to solve Vandermonde systems 

1) Given {j k } compute q={q 0 ,---, 1m)- 

2) Given b and q compute 



for 0 < s < M 
3) Compute p as 



I'M 

" Q'M ' 



_JrP{B 
Git) 



This algorithm is equivalent to the following factorization of the inverse of the 
Vandermonde matrix in terms of a diagonal matrix, its transpose V 1 , and a triangular 
Hankel matrix, 



(5.5) 



This description is a particular case of th< 
[20] or close to Vandermonde matrices [f 
results as (see [20, p. 548]) 



i formulae for Lowne 
, Corollary 2.1 p. 157]. We > 



-Vandermonde 
an state those 
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where the vectors x = (x 1; - ■■, x M Y and y = (y 1 , ■ ■ ■ , y M Y are solutions of 

Vx = (0, • ■ ■ , 1)' and V*y = [7^]^. 

Since 7 r are the roots ol' Q(z) we can lake y = — ( (/ 0 . - - • . y.i/- 1 )* and if li{z) = z M in 
(5.4), then P(z) = 1 and x = (^, • • • , ^y) 4 . 

Remark V.l For Algorithm V.l, we first obtained the eigenvector q corresponding 
1,0 an eigenvalue close to e. Thus, Step 1. of the Vandermonde algorithm is already 
accomplished and Step 2. can be performed using the FFT. Furthermore, the nodes 7* 
belong to the unit circle and via t he unequally spaced fasl Kourier transform we have a 
fast algorithm to obtain the weights. 

Remark V.2 As an example, we use this approach to derive the solution of the Vander- 
monde system with nodes at j r = e' 2 ^^*^, 1 < r < M. In this case, Q(z) = 1 — z M and 
P(z) = z M -'P(z-') = Q{z)B(z) + z M ( : --) = -B(z) + z M '{--■). We conclude P = -B(z) 
and thus 

P{lr) = -7 r M - L fl(7r) 

Pt Q'( 7r ) 




,-i2jrr*/M 



As expected, we obtained the inverse of the discrete Fourier transform matrix. 
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VI Approximation of bandlimited functions 

Let us consider the problem stated in (1.3)-(1.5), that is, given 

u(x) = j'w(T)e icTX dr, (6.1) 

construct the function u(x) in (1.4) such that (1.5) holds. We show in this section that 
the approximation obtained with exponential sums holds in any subinterval of (-1. 1). 
In fact, in Theorem VI. 1, we prove the existence of such an approximation even though 
our proof does not provide a practical method to obtain the nodes and weights in the 
exponential sum. In practice we obtain them using Algorithm V.l. 

We assume that we have access to values of u(x) for x uniformly sampled. We 
select the sampling rate to be at least twice the Nyquist sampling rate for u(x). We have 
observed that this is the minimal rate for which our method works properly. 

In fact, let us discretize u(x) at nodes x k = for < N and pick N such that 
JV > The value of N determines our sampling rate. The resulting values are 

u k = u(x k ) = £w{r)4"*dT. (6.2) 

Defining v = then v < 1/2 and by changing variables t = ut, 

Uk = J a (ty* tk dt for |Jb|<JV (6.3) 

are the trigonometric moments of a new weight a(t) = £iw(£) supported in (— v, v) C 

Now, assume we can approximate J_ 1 cy(t)e'" ty dt by J3^ = , Wje l "°' v for \y\ < N. 
Then, since 

u(x) = j\{T)^ r 'dT = £a(t)4' tl "dt, 

we can appproximate u(x) for \x\ < 1. 

Indeed, we now show that for any d, 0 < d < 1, we can approximate u(x) for 
\x\ < d. 

Theorem VI. 1 Let a be a weight supported in [—v, v], 0 < v < | and let e and d be 
posilivt itunibtrs until d < I. I'litu, for \ mjJkieiiUy laiyt, Ihtrt vjisl nal comlanh 
{wi , ■ • • , wn} and {9\ , • • • , 9 N }, with Wj > 0 and \0j\ < v such that 



(6.4) 
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For the proof, we will use the fact that exponential functions can be well approx- 
imated by bplines interpolating them at integer values. 

For fixed t e [— tt, ty] and positive integer m, consider the exponential Euler spline 
of order 2m - 1, 

^i(:^)=£er«Wi(s-*), (6.5) 
kez 

where L n (x) is the fundamental cardinal spline of order n, L n (r) = 5 ra , for all r£Z. 

We will use the following properties (see [21, pages 29, 30, and 35]) valid for all 

x, t e R, 

\S 2m -i(x,e u )\ < 1, (6.6) 

|e i ^-5 2m _ 1 ( 2 ;,e irt )|<3|i| 2 " 1 , (6.7) 

and <d m er^, (6.8) 

for positive constants d m and a m . 
Proof of Theorem VI. 1 Let 

u(y) = J a(t)J'»dt 
and for each rrc define the spline of order 2m — 1 interpolating u(y) at the integers, 
a(y) = 



^ ff(t)ft m -i(»,e I,rt )A. 

By (6.7), 



!«(»)- «•(»)!< 3 



y*ff(t) |t|*"dt<3i^»|H|i, 



where \\a\\ L = J'_^ cr(i) efc. We choose m such that 3i/ ,im ||<7||i < 

On the other hand, for each N, Theorem III. 7 allow us to represent the moments 

u(k), \k\ < N, 

u(k) = J a{t)e i * kt dt = J2 VJ j >iiwe ' k +"<>(- 1 ) k > ( 6 - 9 ) 

whprp 

4|kl 



" 2 + (2-V3) AI + (2- 



(6.10) 
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fi(y) = X>e>"" 
then u(k) = u{k)+w 0 (-1)* for |A;| < N, and defining 

a(y) = £fi(*)Wi(v-*) 

i=i 

(6.7) gives the estimate 

|%) - S(y)| < 3^^l^| 2m < 3v*"(u(0) - w 0 ) < 3^'Hall! < |. 

We have shown that u(y) is close to a(ij) and u(y) to fi(j/). To finish the proof, we need 
to show that \a(y) - a(y)\ < f , for \y\ < dN + 1. Now, 

<«(») - a(») = £ w ° (-1)* wi (»-*) + E («(*) - «(*)) wi (» - fc) 

|*|<JV |*|>JV 

= Wo S 2m _ 1 ( 3 /,e")+ puttJ-iW-ml-ljVi&f-*) 

and 

|m(A;) - w(fc) -w 0 (-l) fe | < K*0| + |«(AOI+"o<f> + f>+"o 

.7=0 J=l 

< 2u(0)=2||<7|| 1 , 

where we used (6.9). 

Using (6.6) and (6.8), 

Ww)-fi(y)l < + 2 Hi ^£e^-l^*L 

|*|>JV 

Hence, for large N, we can estimate tu 0 using (6.10) and for the last sum, when \y\ < 
dN + 1, 

^ e -a m |»-*| _ ^ ( e -a m (*-») +e -«.W) 
|*|»V 
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Remark VI. 1 In the proof of Theorem VI. 1 we used Theorem III.7 to represent the 
moments (6.3) as 

J a(t) e^ tk dt = J2 Wj e l ' 0 > k + w 0 (-l)*, 

where w 0 dec:reases exponentially with ,V. The approximation in (6.4) is obtained using 
these Wj and 9j. Nevertheless, in practice, we use instead the weights and nodes from the 
quadratures in (4.12). We also note that, as the proof of the theorem indicates, for (1.5) 
to hold on most of the interval (-1, 1), we should appropriately oversample the function 
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VII Approximation of integrals by linear combina- 
tions of exponentials 

In this section we show that by solving the problem (1.3)-(l.o) for x e [—1, 1], we also find 
quadratures for functions that can be considered as a linear combination of exponentials 
with bandlimit c. We provide two examples with Bessel functions and with PSWF. 

Proposition VII. 1 For x £ [—1,1] let us consider- 

V( X ) = j\(T)J 2n (cXT)dT, (7.1) 

where w > 0 is a weight, J' ln is the Bessel function of order In, n > 0 and c is a positive 
real constant. Then we have 

\v(x)-v(x)\<-e, (7.2) 
where ^ 

v(x) = f2w k J 2n (cx6 k ), (7.3) 

and the nodes 0* and the weights Wk are as in (1.4) but for the weight w defined as 
w(t) = w(t)/2 for 0 < t < 1 and w(r) = w(-r)/2, -1 < r < 0. 

Since J- in is an even function, we have 

V{x) = j w(T)J 2n {cXT)dT. (7.4) 

Using 

J 2n (0= { -=±r J^J^jytdy, ( 7 .5) 

we obtain 

v{x) _ Hx) = t2L £ [ 1' w{r)fdT - dy. (7.6) 

Since \y\ < 1, we have (by selecting nodes and weights as in (1.4)) 
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. Thus, we obtain 



\v(x) - v{x)\ < ^ ^M^dy = ( - jT | cos(2na:)|da; = { f 6 % ^ (7-8) 

A similar result holds for the PSWF which are defined as the eigenfunctions of 
the operator 

F e {4>)(x) = (7.9) 

where c is a positive real constant (bandlimit) and F c : U 2 [— 1, 1] -» L' 1 [— 1, 1]. These 
bandlimited functions satisfy 

Wj{ x ) = j e icct ipj(t)dt, (7.10) 



where the eigenvalues Xj, j = 0, 1, . . ., are all non-zero and simple, and are arranged so 
that |Aj_i| > \Xj\, j = l,2,.... 

Proposition VII.2 For- non-negative integers j, let us consider integrals 

1 w(r)v ;] (r)dT, (7.11) 



where w > 0 w a weight and ij> } is the PSWF corresponding to Ifu band Until c > 0. 
T/teri 

K-^|<^ €! (7.12) 



fy^toWiWk), (7.13) 
and the nodes 9^ and the weights are the same as in (1.4). 

For large c, the spectrum of F c can be divided into three groups. The first group 
contains approximately 2c/7r eigenvalues with absolute value very close to one. They 
are followed by order log c eigenvalues whose absolute values make an exponentially fast 
transition from I to 0. The third group consists of exponentially decaying eigenvalues 
that are very close to zero. For precise statements see [23], [24], [13] and [28]. 

Therefore, it follows from (7.12) that, for the first w 2c/w eigenfunctions, the 
integrals in (7.11) are well approximated by the quadratures in (7.13). To prove (7.1.2), 
use (7.10), to write 



(7.14) 
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Sir 



1*1 < 1. 



have 




(7.15) 



and ||V>j|| 2 = 1 implies |i/>j(f)|d* < v^- □ 

Remark VII. 1 If the weight luriclkm w is chosen as in (8.40), then the eigenpolynomials 
are the discrete PSWF [22] (see Example 1 in Section IV.3). In this case, the nodes {9 k } 
are zeros of a discrete PSWF corresponding to a small eigenvalue. Therefore, these nodes 
are Gaussian nodes for PSWF (appropriately scaled to the interval [-1, 1]) as stated in 
Proposition VII.2. 



Gfiimhm quailrnluras for exponentials 



40 



VIII Interpolating Bases for Bandlimited Functions 

In this section we construct bases for bandlimited functions with bandlimit c> 0 using 
exponentials {e"*'*}^,, where {t ( } are some quadrature nodes with \t t \ < 1. In particular, 
we derive interpolating bases as linear combinations of such exponentials. 

We start by constructing, for some e > 0 and bandlimit 2c > 0, the nodes |tf| < 1 
and the weights w t > 0, I = 1, . . . , M, where M = M(c, e), such that for all a; £ [-1, 1], 

|/V^-i>e 2iCTtl |<e 2 , (8.1) 

where, Cor each / i hero exisls /' such lhal I.,, = -/., and ir v = «>,. 
Since e 2icxt dt = sin(2car)/cx, we have from (8.1) 



sinc^-*) " ic(x-t)J < g 



(8.2) 



where |x|, |i| < 1. 

In considering bandlimited functions we will use the PSWF (see [23], [14], and 
a more recent paper [28]). The PSWF are real eigenfunctions of the operator F c in 
(7.9) with eigenvalues \j, j = 0,1,.. ., such that |A 0 | > |Ai| > . . . > 0. They are also 
eigenfunctions of the operator Q c = ^F*F C , namely, 



with eigenvalues 



(8.3) 
(8-4) 



Theorem VIII.l For x in [-1, 1] and for any \b\ < c and e > 0 there exis 
{ a i)iLi an d constants Ay, A 2 such that 

||e ite -^a,e ict|X || <A l( , (8.5) 

and 

ije^-^aie"''*! < A 2 e, (8.6) 

where the nodes \ti\ < 1, I = 1, . . . ,M are the same as in (8.1) and do not depend on b. 
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In other words, for a given precision e, the functions {e"*'*}^! form an approximate basis 
for bandlimited functions. 

The proof of this theorem uses the fact that the PSWF are uniformly bounded 

IMIU^CU J = 0,1,..., (8.7) 

where C x does not depend on j. Through direct numerical examination it is not difficult 
to verify that is a fairly small constant that weakly depends on the bandlimit c. 
However, we are not aware of a proof that provides a tight bound in (8.7). A proof that 
C x exists can be constructed using the fact that for j >- c, the functions il'j approach 
the Legendre polynomials. In order to obtain appropriate estimates, one can use the 
recurrence relations derived in [28]. 

Proof of Theorem VIII. 1 Let us start by expanding e ita into the basis {vj\f',j 
corresponding to the bandlimit c. We have 

3=0 

where \b/c\ < 1, and by (8.7) 

U _ X)A^(6/c)^(x) < jt \M b / c )\ I^WI < C x J2 |A,| \Mz)\. (8.9) 

j=0 I j=M j=M 

Thus we obtain 

||e ita -f^W^(z)|| <CS,EW, (8.10) 
and, using (8.8) and the orthonormality of fj on [—1. 1], 




From (7.10), we have 

J e- ictt 'il> j (t)dt = \ j rp j (t,), (8.12) 
and, using (8.2)-(8.3) and the relationship between and A.,- in (8.4), we obtain 

M ! 

g^e^A^t,) - |A/^(z) < e 2 J i \^(t)\dt < %/2e 2 . (8.13) 
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Thus, we arrive a 



£wie^^j(ti)-^(i) <2^. (8.15) 
Combining (8.10) and (8.14), we have 

_ £ £ Wie -V 3 -(6/c)^(ii)| M + v^C £ ^7y ( 8 - 16 ) 
Similarly, combining (8.11) and (8.15), we obtain 



M M-l II 00 M-L 2 

"EE ^e^(Wift) < Ccc J £ |A,P + 2 Coo £ ^, (8.17) 
By setting 

a, = io,^ ^(&/c)^(t,), (8.18) 
j=o 

and observing that |A M | « « and that |A,| < |A M | for j > M , we obtain (8.5) and (8.6). 
□ 

We now construct two useful bases as linear combinations of the functions {e ict ' x }fl , . 
First, let us consider the following algebraic eigenvalue problem: 

f2^ ct "' t ^ j (t,)=V^j(tm), (8-19) 
(=i 

where t ( and w t are the same as in (8.1). By solving (8.19), we find r)j and We 
then consider functions j = 1, . . . . M, defined for any x as 

1 M 

ttyx) = -^tu^^fo). (8.20) 
^ (=i 

'Hie functions in (8.20) are linear combinations of the exponentials {e i<:x '< }" =1 . We 
will show that the functions *j are nearly orthonormal and we will use them as an 
approximate basis for weighted bandlimited functions with bandlimit c. 
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Remark VIII. 1 The functions *j mimic the PSWF. However, one has to be careful 
relating -ipj and 9j. Since a large number of eigenvalues A ; satisfy both Xj is close to 
\/2ir/c and rjj - Xj = 0(e 2 ), then, in solving (8.19), we also obtain a large number of 
eigenvalues r\j with absolute value close to y / 27r/c. Therefore, in spite of the fact that 
all Xj are distinct, we may obtain a group of eigenvalues that are identical within the 
precision of the computation. In this case the functions <l' y correspond to linear combi- 
nations of the PSWF. In other words, we need to impose additional conditions in (8.19) 
to maintain a proper correspondence with the PSWF. However, in many applications 
there is no apparent need to make such an identilication since in all cases the resulting 
functions span the same subspace. 

Proposition VIII. 1 The functions tyj, j = 1, . . . , M, are nearly orthogonal and satisfy 



Proof We start by denning qj = ^/wi^fj(ti). We substitute in (8.19) to obtain 



The vectors {<f J } are eigenvectors of the matrix S, 5 m ; = yw^e' 1 ' 4 '"*' y/wi. If we take into 
account the symmetry of nodes ti and weights wi in (8.1), we obtain that the matrix S 
is normal, SS* = S*S, and, in addition, S = S*. In Proposition VIII. 2 we show that for 
such matrices there exists an orthonormal basis of real eigenvectors. Thus, computed via 
(8.22), we assume qj to be a real orthogonal matrix and then 




(8.21) 



M 



Vw^e wtmt 'V^Qi = m qL- 



(8.22) 



(8.23) 



and 




(8.24) 



We have 




dt 



(8.25) 
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and, from (8.1), we obtain 




M M 



(8.26) 



For the last inequality we also used Schwarz inequality and (8.24) to estimate 



To finish, using (8.19) and (8.24) we simplify (8.26) and arrive at (8.21). □ 
We still need to prove 

Proposition VIII. 2 Let S be a normal matrix such that S = S*. Then there exist an 
orthonormal basis of real eigenvectors of S. 

Proof First, since S is normal, eigenspa.ces corresponding lo dilTcrcnl eigenvalues are 
orthogonal. Thus, it is sufficient to prove the proposition for any eigenspaee /-.'( A) = 
{* : Sx = \x}. Also, normality of ,9 implies that for any x € E(X) we have * £ E(X). 
Indeed, from S*x = Xx and S = S* it follows that Sx = Xx. Consequently, if {v k }™ =1 
is a basis of E(X), then E(X) can be spanned by the real and imaginary parts of v k . 



where for any vector x = (ii, . . . , xm) : Re(x) — (Re(x,), . . . , Re(x TO )), and similarly for 
Im(a:). By Gram-Schmidt orthonormalization of the 2m (linearly dependent) vectors A, 
we obtain the desired result. See another proof in [9, Theorem 4.4.7] □ 

Let us now construct interpolating bases as linear combinations of the exponentials 
{e icxt '}f =v We define functions B. k , k = 1. . . . , M as 



!>,«>,< ^(to*^) < E u "i*^)iE u "i*j'(*»)i 

< (f> |^(*,)|) 2 < (f>)(f> = £>. 



A={Re(v k )Am(v k )}f =l 



(8.27) 



where 



m 1 M 1 



(8.28) 
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By direct evaluation in (8.19) and (8.23), we verify that functions R k are interpolating, 
Rk(t m ) = 5 km . (8.29) 
Let us show that the integration of R k (t) e' at , where \a\ < c, yields a one point quadrature 
rule of accuracy 0(e). 
Proposition VIII.3 For \a\ < c let 

A fc = y R k (t)e iat dt-w k e iat K (8.30) 

Then we have 

where ||A|| 2 = ^Eili |A*| a . 

Proof For functions e mx , where \a\ < c, we have 

jf' J4(i)e^ di - nu f; ii^er*-*^ = f> iSj , (8.32) 

where 

,i m 

Sl = j ^ t >+^)dt-Y / ^ ctm{t,+a/c) - ( 8 - 33 ) 

We obtain using (8.27) and (8.29), 

E^E;^ 1 = £ ^^(^)e-- = (8.34) 

and, therefore, A k in (8.30) can be written as a matrix-vector multiplication, A k = 
Hiii r kii>i- The inequality (8.31) is then obtained via the usual Z 2 -nomi estimates, taking 
into account that the matrices qj, and qj in (8.28) are orthogonal. □ 

We have observed (via computation) that max* = i,...,M \w k \ = O(l) and mhn—i,...^ |%| = 
0(e _1 ) in (8.31), thus resulting in || A|| 2 = 0(e). Next we derive a weak estimate showing 
that the functions R k are close to be an interpolating basis for bandlimited exponentials. 
Proposition VIII. 4 For every b, \b\ < c let us consider the function 

n b (t)=e iM -jrj Uk RM (8-35) 

k=l 

Then, for every \a\ < c, we have 

f n.(t) dt\ < (1 + M m ^='."-,«M ) e*. (8.36) 
|7-i | min* =1 ,.., M |%| 
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Proof Using (8.30), we have 

J fi»(t) e iaf dt = J e^' dt-Y^w k e'<*** - ^ e iM * A fc , (8.37) 

where 

A k = j B. k (t)e iat dt - w k e iat " . (8.38) 

Applying (8.1), we obtain 

n»(t)ef**| < e 2 + %/M||A|| 2 . (8.39) 
The estimate (8.36) then follows from Proposition VIII.3. □ 

Remark VIII.2 Using the functions l! k , k = I..... M, on a hierarchy of intervals, it 
is possible to construct a multiresolution basis (for a finite number of scales) similar to 
multiwavelet bases. We will consider such construction and its applications elsewhere. 



VIII. 1 Examples 

For the weight 

"<*>={! ll!*Z ] : a - k ^ 

we construct a 30-node quadrature formula so that (8.1) is satisfied with e 2 ss 10~ 15 . 
We display the error in Figure 9. For bandlimit c = 7.5 n, we construct an interpolating 
basis {if*}!!, and display one of the functions in Figure 10. We then demonstrate the 
performance of the interpolation using this basis for three examples, 

COB(d), (8.41) 

gait) = t, (8.42) 
flj(t) = P 9 (t) , (8.43) 

where P 9 is the Legendre polynomial of degree 9. These three functions are not periodic 

and we use 

g 1 (t) = f^cos(ct i )i?,W: (8.44) 

30 

h{t) = 52t,R,(t), and (8.45) 
(=i 

30 

~g 3 (t) = 1 £P»(t l )Ri(t) (8.46) 
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as approximations. We display the function r/i in Figure 11 and the error of approxi- 
mation by gi in Figure 12. Similarly, we display the function <j 2 in Figure 13 and the 
error of approximation by g- 2 in Figure 14, the function g 3 in Figure 15 and the error of 
approximation by g- s in Figure 16. 
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Figure 11: Function /^(t) on the interval [—1,1]- 
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Figure 15: Function g 3 (t) on the interval [—1,1]. 
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Figure 16: Difference g 3 (t) - g 3 (t) on the interval [—1,1] 
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IX Conclusions 

In this paper we have introduced a new family of Gaussian-type quadratures for weighted 
integrals of exponentials. These quadratures are parameterized by the eigenvalues of 
the positive definite Toeplitz matrix constructed from the trigonometric moments of 
the weight. The eigenvalues of this matrix accumulate towards zero and appear to have 
exponential decay (see Figure 1). For small eigenvalues, these quadratures are of practical 
interest. 

The remarkable foal lire of these quadrat ures is thai 1 hoy have nodes outside I he 
support of the measure and, as il. turns out . t he corresponding weights are negative ami 
small, roughly of the size of 1 lie eigenvalue. The case corresponding to the smallest 
eigenvalue is equivalent to the classical Caratheodory representation. 

As an application of the new quadratures, we show how to approximate and inte- 
grate several (essentially) band li mil ed functions. We also have constructed using quadra- 
ture nodes and for a given precision, an interpolating basis for bandlimited functions on 
an interval. 

In the paper we made a number of observations for which we do not have proofs. 
Let us finish by stating two unresolved issues. First, it is desirable to have tight uniform 
estimates for the L^— norm of the PSWF (with a fixed bandlimiting constant) or, ideally, 
for the eigenfunctions associated with more general weights. Second, we conjecture that 
in Theorem IV. 1, it is not necessary to require distinct roots for the eigenpolynomial 
since it might be a consequence of the eigenvalue being simple. We neither have a proof 
nor a counterexample at this time. 
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A Appendix 

To prove Theorem II. 2, we use a technique that goes back to [2] (see [27, Theorem 7.3] 
and [18, Chapter 5] for more details) which involves the Fejer kernel, 

for real x. 

We need the following result. 

Theorem A.l ([18], Theorem 8, Chapter 5) For \k\ < N let 

M 
3=1 

where Pj > 0 and \ Zj \ = 1. Then, for all L, 0 < L < N, 

(L + l)||p||l<cg + 2£>| 2 . 

k=\ 

Proof Let a* = 1 - j)^- be the coefficients of the Fejer kernel Fj, and write Zj = e m9j . 
Since pj > 0 and F L (6) > 0 for all 6, 

5^ dk\ck\ 2 = ^ a k Y^PjPi(f) k 

\k\<L \k\<L j,l 1 

M 

= Y,m F dOj-Oi)>FiXV)Y,p 2 i 

3,1 3=1 

= (£+i)XX 

The theorem follows because a 0 = 1 and a k < 1. □ 

Proof of Theorem II. 2 We first use (2.1) to extend the definition of c k as C- k = 
for A; = 1, • • • , N and c 0 = ft- We then define the Toeplitz matrix T N , (T N ) kj = 
i c :i-k)o<k,i<N anl ^ tne polynomial 



Q(,)=n^-e^) = E^- 
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Then, q = (q 0 , ■ • • , q M , 0, • • • , 0)' belongs to the null subspace of T N because for all / 

M M M 

since, for all j, Q(e m6 i) = 0. 

The matrix A = T N - cqI has the eigenvalue -c u because T N is singular. We 
can bound Cq by the Frobenius norm of A to obtain 

c§<2^i| CA , +w | 2 <2 J V||c||2. 

To finish the proof, we use Theorem A.l with L = N. □ 
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Abstract 

Polynomials arc one ol* the principal loop, of classical numerical analysis. 
When a function needs lo be interpolated. integrated, dillerenlialed. etc, it is 
assumed to be approximated by a p.,h riotnial ol' a certain hxed order though the 
polynomial is almost never constructed expheilly ;. am. I a i.rcaliiieiit appropriate 
to. such a polynomial is applied. We introduce analogous techniques based 
on the assumption that the function to be dealt with is band-limited, and use 
the well developed apparatus of prolate spheroidal wavefunctions to construct 
quadratures, t at cqs <lal i on aad dlflerelillatlon Iota I HI lac, etc. tor banddimited 
functions. Since haiiddiinued lunclimis are m'lcn eiieoumered in physics, 
engineering, statistics, etc, lite apparatus we introduce appears to be namral 
in many environments. Our results are illustrated with several numerical 
examples. 



1. Introduction 

Numerical quadrature and interpolation are a well developed pari oi numerical analysis; 
polynomials are the classical tool for the design of such schemes. Conceptually speaking, 
one assumes that the function is well approximated by expressions of the form 

7=0 

with reasonably small is, and one deigns algorithms thai are effective for functions of the 
form ( I ) needless lo say, one almost never actually computes the coefficients one only 
uses the fact ol" their existence. Obviously, the polynomial approach is only effective for 

functions ihal are well npproxi led by polynomials. 

When one has to handle function-, that are well heha\cd on the whole line (for example, 
m signal processmgl, polynomials are nol an appropriate tool. In such eases, trigonometric 
polynomials are used: existing look are very satisfactory lor dealing with fmietioiis debited 



and well behaved on the whole of E 1 . such tools, in effect, make the assumption that the 
functions arc band-limited or nearly so: a function / : R -+ K is said to be band-limited if 
thereexisl i isitive real id a function cr e L -[- 1 . 1 ] such that 



However, in many cases, we are confronted with band-limited functions defined on intervals 
(or, more generally, on ci impact rogn ns in Wave phenomena arc a rich source of such 
function bolliiiitheengiiiecnii i i uta wK ntexl Ills iu Iso encountered in lluid 
dynamics. signal processing, and many other areas. Often, such functions can be effectively 
approximated by polynomials via standard Mils of classical analysis. However, even when 
such approximations arc feasible, they arc usually not optimal. Smooth periodic functions 
are a good illustration ol this observation: while lhc> run be approximated b\ polynomials 
(for example, via ( Tiebyshev or Legendre expansions), they are more efficiently approximated 
by Fourier expansions, belli h a" analytical ;uid numerical purposes. It would appear that an 
approach explicitly based on trigonometric polynomials could be more efficient in dealing with 
band-limited functions. 

In the engineering context, such an apparatus was constructed more than all years ago 
(see |8 1(1,21,221). The natural tool for ana lysi ng hand-limned functions on R 1 is the Fourier 
trans ion n, unless the 1 unci ions arc periodic, in which case 1 he natural tool is the Fourier scries. 
The authors of |2 1,221 observe that, for the analysis of band-limited functions on the interval, 
profile spheroidal wa vol unci ions . I 'SWT s; arc likewise a natural approach. I he authors also 
construct a nuilli-dimensioriul version of the theory. Touch their apparatus is only complete 
for the case of spherical regions. 

The present paper constructs tools lor using the approach of [21,22] in the modern 
computational environment. We construct a class oT quadratures fa band-limited functions 
that closely parallel the tlaussian quadratures lor polynomials. The nodes are very close- 
to being roots of appropriately chosen PSWFs. the resulting quadratures arc stable, and all 
weights are positive. As in the case of polynomials, there arc interpolation. Ji flerenliation and 
indefinite integration schemes associated v. it h l he ohtai ned quadratures, exact on certai n classes 
of hand-limited functions. These procedures are the main tools necessary for the numerical 
use ol' spectral discretizations based on PSWFs. instead of on the usual polynomial bases. 
When dealing with band-Inn hed functions, the number of nodes required by these procedures 
to obtain a prescribed accuracy is much less than that required by their polynomial-based 
coiuncrparls. An additional bonus is the fad thai the condition number of differentiation 
of PSWFs is less than that ol dillereiilialion of the usual polynomial basis functions (see 



This paper is organized as Pillows. Section 2 suiuinari/.es various standard mathematical 
facts used in the remainder of the paper. Section a contains derivations of various results used 
in the algorithms described in later sections. Section 4 describes algorithms for evaluation 
of PSWFs and associated eigenvalues. Section 5 describes a construction of quadratures 
for band-limited functions. Section 0 describes an alternative apiproach to arriving at such 
quadratures: it shows that roots of appropriately chosen PSWFs can serve as quadrature 
nodes. Section 7 analyses the use of I'SWFs lor interpolation. Section 8 contains results 
of our numerical experiments with quadratures and interpolation. Section 9 contains a 
number of miscellaneous properties of PSWFs, and section 10 contains generalizations and 




(2) 



.8). 



2. Mathematical preliminaries 



this paper the norm of a function is. unless stated otherwise, il 
\f(x)\ 2 dx. (3 



2. 1. Chebyslm- systems 

IKIinilion 2.1. ,1 my/ho/, e offiou lions cJ| <J>„ will he referral to as a Chebyshev sy 

on the interval |«. /-I il emit oj them is ennlinuous ami the determinant 



|«/>„fe) ... <ji„(x„) I 

is noieetofor tinx setjuenee ofpoinl\ .\ \„ su< h thai a " i, --. .i^--. --. r„ ■: /<. 

An alternative definition of a Chebyshev system is thai any linear combination of the functions 

Hxarnples of Chebyshev and extended Chebyshev systems include the following. 
(Additional examples can be found in [12].) 

Example 2.1. The powers 1, x, jr 2 x" Jorm an extended Chebyshev system on the interval 

(-OO. 00). 

Example 2.2. ihe exponentials e ' 1 ' . e e ' /roan mi extended Chebyshev system 

for any Xi X„ > 0 on r/it- interval [0, oo). 

Example 2.3. TTie functions- I . cos x, sin x, cos 2i, sin 2x cos nx, sin nx form a Cheby- 
shev system on the interval [0, 2ir], 

2.2. Generalized Gaussian quadratures 

A quadrature rule is an expression of the form 

5>7*C»J>. (5) 

where the points Xj eM. and coefficients Wj e R are referred to as the nodes and weights of 
the quadrature, respective!}. Ibex serve as approximations to integrals of the form 



with m being an irileyrable non-negative function. 

Quadratures are lyjncally chosen so thai the quadrature 10) is equal to the desired 
integral (6) for some set of functions, cvmmi.nl> polynomials of some fixed order. Of these, 
the classical Gaussian quadrature rules consist of n nodes and integrate polynomials of order 
In - 1 exactly. In [14], the notion of a Gaussian quadrature was generalized as follows. 
Definition 2.2. A quadrature formula will he referred to as Gaussian with respect to a set of 
2n functions </>,,..., <h„ : [a, b] -* R and a weight function at : [a, b] -> R* if it consists 
i / i, Mi 0 a , 1 , ■ h aw O'o "os - ' , o ., a, , , o,. iion", a'' I ,V /< u, o, ., , / - 

.a// i 1 2«. ///(' weights and node', oj a Gaussum (luudialuex veil I V iclcned lo a\ 

Gaussian \yei'.;hls and nodes, resne, lively. 



The following tl n I Mil [15.16] pioofsofitcanalsobe found 

in [13] and [12] (in a somewhat different form). 



Theorem 2.1. Suppose that the functions <tn , . . . , <fe, : \a, h | -> Mform a Chebyshev system 
on [a,b]. Suppose in addition that m : [a,b] -* R is a non-negative integrable function 
[a, b] -> R. Then there exists a unique Gaussian quadrature for the functions fa , . . . , 4> 2 „ on 
\a. h\ with respect to the weight function The »v/«/lO of this quadrature are positive. 

Wh i le the ex i sicnce. .f genera I i/od Gauss in n q i indral i ircsw as observed more than 100 years 
ago, the eonstruelions found in |7. 11-13, 15, 16|. do not easily yield numerical algorithms 

ill il I > i i I t ueled reeenll 

(see [2, 14,27]). 

Remark 2.1. It might be worthwhile to observe here thai when a generalized Gaussian 
quadrature is to be constructed, the dclcrrninaliiai uf its males tends to be the critical step 
illi li lie ] i c lii I | 14 | 1c ii u lie i i i It neou ly). Indeed, 

once the nodes .v,. ay v„ ha\c been lound. the weights uy. uy uy can be determined 

easily as llic solution of the „ /. it system ol linear equations 

Y i wj-Mx j ) = j co(x)Mx)<ix, (7) 



2.3. 1 egeiuhe polynomial-. 



with the initial conditions 

P 0 W=I, Pi(x)=x; ( 
as is well known, 

ft(l) = l (1 

for all A- t), 1,2 and each of the polynomials /', satisfies the differential equation 

„ 2 .d 2 ftM „ dft(.v) 



l)f\ (*) = (). 



£ (P „ ( , ))2d , = _L 



the normalized version of the Legendre polynomials will be denoted by P„. so that 
T n {x) = Pn(x)-Jn+\/2. 



Lemma 2.2. For all integer k > n, 

f • K;, 

For all integer 0 s$ k < n, 

|^ AT* T^W ibc^O. 

2.4. Comolutlonal Volicrra equations 
A omvolutioiml Voltcrni equation of the wcond kind is an expression of the form 

<p(x) = j* K(x-t)<p(t)dl + <,(x) (16) 

squarc-intcgrablc and (/; : [<;. /;] C is the function to be determined. IVoofs of the following 
theorem can ho fumd in 4|. as well as m many hUilt sources. 

Theorem 2.3. liquation tl6i always has a unique solution on tin: interval [a. !>]. If both 
functions K.n are k limes continuously differculialilc. lite solution e; is also k limes 
continuously differentiate. 



(14) 
(15) 



In [his subsection, wo summarize certain lads aKn.il [ho PSWhs. I 'nless slated otherwise, all 
these facts can be found in [18,21]. 

(.liven a real ( v 0, we will denote by /; the operator f. : | 1.11 • /. : | 1. 1 I delmed by 
the formula 



>r all natural ;'. 
esponding to a 



Obviously, F r is c.iinp.ici; we will denote by >. r , , X /.„, . . . the eigenvalues of F r ordered 

so that > \X j \ for all natural ;'. For each non-negative integer j, we will denote by i/fj 
lite cigenlunclions c.»rrcN|'»mlii)g In k ,. so that 



Theorem 2.4. For any positive real ,. the civcnfunciions oo. f l , of the operator F c 

eigenftmctions are even, and the ocid-niiinhcred ones are odd. Alt eigenvalues of F c are 
non-zero and simple: tin- nrn-nmnbeivti ca'cnvtbuco are purely real ami the odd-numbered 
ones arc pin ch n i I i ( / i 

system on tin- interval | - 1 . I ]; in pa ill: a Ian l he Inn, lion <t h,i\ c\uciivi .am. on line unci nil. 
for any i = 0,1 

ol hmaion- on I I . I |. and in lli:u ^.m,"-a. Ihe oownmo [li.n ill. In no ion:. | ■ ■ \ haw nun mem .ni 



We will define the self-adjoint operator Q c : Lr[-1. 1] -+ L 2 [-L 1] by the formula 

u*=iL mc *-7 tyi «™ (19) 

that Q r has the same dg«iilimv-li>iii» us I ... and llial (lie ./Hi cm descending order) eigenvalue 
Hi of Gc is connected w ith , by the formula 

;'i = ^-|>-y| 2 - (21) 

The operator Q c is obviously closely related to the operator P c : L 2 [-oc, oo] [-co, oo] 
dclincil by the formula 

^f<P) = ^ ■ / sm(r ^~ n) . ^(n dr. (22) 

which, as is v.ell known, is Ihc ..nhoponul projection operator onto the space of functions of 

I -'or large . . Ilie spectrum .f n consists of three parts: about 2c ,',T eigenvalues that are 
very close to I, followed by order logo > eigenvalues which decay exponentially from I to 
nearly 0; the remaining eigenvalues are ah very close to zero. The following Iheorem, proven 
i in a slightly different form i in |2'*i|, describes the spectrum ,.f Q. more precisely. 

Theorem 2.5. For any positive real c and 0 < a < 1 the number N of eigenvalues of the 
operator Q c thai are greater than a satisfies the equation 

N = ^- + (J^ log j log(c) + Oflog(c)). (23) 

By a remarkable coincidence, the elgelililnclioiis ,. o i/o. of the operator ij Una 

■ Hit to be the PSWhs. well known In™ classical maihenialieal physics 1 see. e.g.. [171!. Ihe 

following Iheorem formali/es ilus statement: h is prove a eonsiderabl v more ceneral form 

in|22|. 

Theorem 2.6. For any . > 0, tlicre exists a strictly imreasiiix sequent c of [Hisitire real 
numbers xo- Xi, ■ ■ ■ such that,, for each j > 0, the differential equation 

(1 - x 2 ) f"(x) - 2.x f(x) + ' Xj - c 2 x 2 ) ir(x) = 0 (24) 
IMS " solution that is continuous till the interval | 1.1]. l-nr ea. I, ; i ., Ihe function ./, , 

(defined in theorem 2.4) is the solution of (24). 



3. Analytical apparatus 

3.1. Prolate series 

Since the functions V'n, V'f <Pn. ■ ■ ■ are a complete orthonormal basis in L 2 {— 1, 1], any 

formula for the inner product of PSWFs with another function / is also a formula for the 
coefficients of an expansion of /' into PSWhs i. w inch we will refer to as Ihe prolate expansion 
of /). Thus the following theorem provides ihe c, >eflicienis of the prolate expansion of the 
derivative of a prolate spheroidal finicuon. and also ihe coefficients of the prolate expansion of 
a PSWP 'multiplied by .v. '1 hose eoeflicienls are also lite entries ol the matrix for differentiation 



f a prolate exp ai u ci 1 I pa i u 1 the entries of the matrix for 

multiplication of a prolate expansion by i , rcspec lively. . I tiese formulae are not, however, 

utable for prod u i i ^ ses they exhibit c t sti pliic 
cancellation.) 



J 1 f'„(.x)f m (.x)dx = J* > 



11 and odd. id?) is obvious. In order lo 



K 4' n = j e krt t/r„(Od« 
•cliim 2-? I- Diffcrcntiatinp ( 2 8 p with respect to a . wc obtain 
K f,; (x) = i c J t e iri ' (t) dt. 
h sides of (29) on i/<,„ and using 

K, j | \l>'„ (.x) V/„, (x) dx = ^f t V'™ t* )ft c 1 "' f„ (<) oV d* 
= i £ j | / V'„ (0 f ^ e iOT V m (*) At di 
= icl m j ( lA,»A.(l)dl. (30) 
Obviously, the above calculation can be repeated vvilh «; and a exchanged, vielding the identity 
A™ | ^ *'„, (x) f„(x) dx = ic>.nf i < iMO (0 df ; (31) 
combining (30) with (3 1 ), wc have 

J ^ M ifc. « dx = jff^™ M +'» W d-T. (32) 
On the other hand, integrating the left side of (32) by parts, we have 
j i K to Va, (x) dx = V/,„ ( I ) ifc, (l) - <//,» i-DW-D- / 1 € to i'm (x) dx. (33) 



Sincem ^ n(mod2), v 



/ i v/,;, (x) v/„ w it = 2 v/„ ( i ) f„ ( i) - y f n o 



Now, combining (32 j and (34) and rearranging terms, we get 

j i fc(*)d* = ^f^lMl)^(l)- (35) 

y * yr„ (*) ^ ) dx = Y c y- j t'„ « f m (*) dx 

= | I ^7^(l)V^(l)- (36) 
□ 

'I Ik following corollary, which is an immediate conseij nencc of i 32 i, finds use in the 
numerical evaluation of the eigen\ allies [a.). 

Corollary 3.2. Suppose lint! c is real tiini positive, ami that liw iitwgtrs m oriel n are non- 
negative. Ifm^tn (mod 2), then 

^ = { l _,*'„(x)* m {x)dx (37) 

X « f^i^(x)t„(x)dx' 

3.2. Ih-rax III i I I I ; / i / 



Lemma 3.3. Let I \ I i i i I i i i / 

£ C '<" T„(x) dx = g £ .v 2 * T«(x) dx + i g & £ * 2t+1 ft(x) dx (39) 



luillierinorejorall integer m L« • MJ + 1. 

| J ' e<" ft « dx - J2 <** / ' -r 2 " ^ to d * "iD« £ -v 2 *'*' ^ to | - 



In particular, if 

n3>2(|e-|«|J + l), (44) 
e'™T„{x)&x\ < (I)"" 1 . (45) 

Proof. Formula (39 i follows immediately from lemma 2.2 and Taylor's expansion of e : ' r . In 
order lo ]irove <4.<). we assume dial, rn is an integer sued l.hal 

m>\e-\a\\ + l. (46) 
Introducing llle notallini 



(47) 



we immediately observe thai, due lo lemma 2.2 and the triangle inequality. 

i^i<£(if-vSr)<£7f- (48) 

Since (46) implies that 

lor all inlegcr m , A . li. we rewrite <4S) as 

mill ohlnm (43! immediately using Stirling's lurunila. lunally, we obtain (4?) by choosing 

in=Le-|fl|J + l- (5D 

□ 

Theorem 3.4. / i ' 

the till nonmilk-rtl l.e.geinlic polynrvmul Ul<;lin<;l in ' /..')/ <«></ /O >..„ /-'<' t/n? eigenvalue which 
mrresimuls lo <fr m lx I («.v fn theorem 2.4 1. J hen lev all iiilcucr m I 0 iw</ o// real positive t, 
ff 

*»2(|*-c| + l), (52) 
I f' fcW AWdjcl < 7— • ( t\ ■ (53) 



Moreover, given any e 



k > 2 (L* • cj + 1 ) + log, Q j + log, , 



H : i'l/></>«-H 



and ruiiicmhcriiig thai llie prolate functions have unit in 
J |V'„W|d.»: < V2, 



3 implies that 



i/>»H«n-Gr/>«'*<ciG: 

Substituting (54) into (53), we immediately see (55). 



4. Numerical evaluation of prolate spheroidal wavefunctions 

Both the classical Bouwkamp algorithm (see. e.g.. |]]) lor the evaluation of the functions (fry, 
and [he ulporilhm preseriled in litis paper for the same (ask. are based on [he expression of 
those functions as a Legendre series of the form 



C* + .*!i2* li 
' ■ at-2 = 0. 



(k + 2)(k + l) 2 , / 2k(k + 1) - 1 2 _ \ , 

(2* + 3)V(2* + 5)(2*+l) ' ' ' P,t+2 V (2* + 3)(2i - 1) ' C *7 ' * 

for [he c. leflieienls /i ; . ;>•• . ... of ihc oxparisn >n 

lor each y II. 1.2 we will denote by p> llie \ eel or 11 /- deli lied by the formula 

f)i = (){,£,...). 

I he following llieorem reslales [he reeursion i 62 •in n slightly different form. 



Theorem 4.1. The . 

formulae 



tries of the matrix being . e 



of the elements of A with even-numbered mis awl columns and Ihe sivnJ consisting of the 
elements ol A Willi odd-numbcrcd rows ami columns. While these two matrices are mtinite, 
antl their ernries d. > n. >l decay much with increasing row or column number, the eigenvectors 
!/;■'] of interest (those corresponding W ihc tirst m profile spheroidal functions! lie almost 
entirely in the leading rows and columns of the matrices (as shown by theorem 3.4). Thus the 
evaluation of prolalc sphenoidal luacllons can be performed by I he follow nig procedure: 

• Generate Ihe leading /-. rows and columns "1 A. where A is given by t54i. 

• Splil Ihe gcncrnled portion of A inlo A^. and .1 , ,,. and use a solver lor the symmetric 
Indiagonal dgcnprnblciii 'such as that in I.APAt Kj t.» compute their cigem eclors {p'\ 
and eigenvalues {/,-}. 

• Use Ihe obtained values ol the coeflicienls (if. ft' . ... hi the expansion l 63 ) to evaluate 
Ihe Ihnelioii ig. al arbitrary points on Ihc inlcrval [-1. I]. 

As a numerical diagouah/aliou of a positive dclmii.e indiagonal mains with well-separated 
eigenvalues, this precoiiipillation slage is numerically robusl and efficienl. requiring Od'w) 
operations to construct the l.egendre expansions oi ihe f nii i 64 1 lor Ihe firsl m prolale 
spheroidal functions; each suhseqiient evaluation of a prolalc spheroidal function takes O(c) 



4.1. Numerical evaluation of eigenvalues 
Although the above algorithm for the eval 
of the differential operator 1 24 1. it does not 
F c (denned in (17)). Some of those eigen 



M-i(v) 



evaluaiing die iniegral on ihe righidiand side numerically; however, that evaluation obviously 
has a condition number of about 1 /a... and is ihus mappropnalc lor computing small kj. A 
well-conditioned procedure is as follows; 
• I 'so (60 1 lo ca leu laic evaluaiing die righidiand side numerically, and with x = 0 (so 
that ijr a (x 1 is not small). 



• (,'se the calculated together with corollan 3.2. to compute the absolute values 
for; 1.2 .. »/. computing each k, from k, 1 1 i and again, evaluating the required 
integral-, numerically). 

• I. 'so the 1'acl that a, i ' a, | i«e theorem 2.4) to finish die compulation. 



Since the PSWl's <//», i/o i», constitute a complete orthonormal basis in L 2 \-l, 1| 

rsec theorem 2.4), 

ei "' = fl ( / L eir " *J (V) dT ) * w< (70) 

for all x, 1 e [-1, 1]; substituting (18) into (70) yields 

<**' = f> (71) 
'linn iffi quadrature inlcgralcs cxnclly the [irM „ cigcnfimctions, thai is, it' 

X>iVcita) = j\j(x)&x.. (72) 

for all j 0, I ii 1. then the error of the quadrature when applied to a function 

|>e-*-/_W 

= X>(£ - £ (f>J *y(x))ir. (73) 



Prom (74), il is obvious that the error of integration 1 73 ! is of roughly the same magnitude as 
A,,, provided that » is in the range where the eigenvalues {A.;} arc decreasing exponentially (as 
is the case for quadratures ol am useful accuracy, sec theorem 2.?) and provided in addition 
Ihal die weights {i/s j are nul large. 

Now, the existence of an n/2-poiiit quadrature that is exact for the first n PSWFs follows 
from the combination of theorems 2. 1 and 2.4; an algorithm for the numerical evaluation of 
nudes and weights of sue]) quadratures can he found in |2|. An alternative procedure for the 
construction of quadrature formulae for band-limited functions (leading to slightly different 
nodes and weights) is described in llie loll owing seel car: a numerical comparison of the two 
can he found in section S Mow. 

Remark 5.1. The above text considers onl 



g(x> = L 1 



||G|| = ||«||(_»..^. ; lhal is 1" suy. allliougli llic arror M such a quadrature when applied lo a 

band-limited function i- n i I" n i led pn ij I fill f that function on the interval 

(il integration, u is h.Hinded proportional to die norm "1 dial function on the entire real line. 

6. Quadrature nodes from roots of prolate functions 

An alternative to the approach of the previous section is to use roots of appropriate PSWFs as 
quadrature nodes, with the weights determined via the procedure described in remark 2,1. The 
following theorems 1 provide a basis fa- hi is; uiuiicricalh i see section X i die resulting quadrature 
nodes tend to be inferior to those produced h\ the optimization scheme of |2, 14,27]; however, 
Ihcy arc useful as starling prints lor dial scheme, or as somewhat less efficient nodes which 
can be computed much more quickly. 

Hie following Kvo theorems conslitute a slraiglillorvvard extension to handd ed functions 

of Kuclid's division! algorithm f.-r polynomials, d'heir proofs are quite simple, and are provided 
here for completeness, since lite author failed lo laid ilium la die literature. 
Theorem 6.1. Suppose that a, <p : [0, 1] C are a pair of c 1 -functions such that 
<PW # 0. 

t is a iH'siiivr veal mnvhev anJ (he f;tm Hens I', ;> are Jelmed <<y llie f, -nualae 



Then titers exist two <•' -functions n. f : [0. 1] ->• C such that 

/W = pW?W + rW 
for all x € R, with the functions q , r : [0, 1] > R defined by the formulae 

r(x) = J^(ty"' dt. 
Proof. Obviously, for any limclions p. ,/ given by . to . and (SI l. 
p(x)q(x) = £ <p(t)e c " dt- J* »,(r)e ic " dr = jf' jf' <p(t) n (t) e^' + '> dt dt 
Defining the new independent variable u by the formula 



I. I hi- -pin 'if inn-ji.irion 



p(x)g(x) = ^ e kM £ <p(u-x )))(!■) >iz du + j e' c " ^<p(u - r) r,(j) dr du (85) 
(sec figure I). Substiiiiting (78i. (82) and I 85 ) into (80 i. we got 

jT e k "^ ¥ >(«-T) ) ;(T)drd« + ^ e™ jf ^(« - r) r,(x) dr du + ?(r)e fa 'dr 

= jT 1 ' 2 fr(Oc 2i "'d< + ^ 1 fT(nc 2k "dr. (86) 



"I Ihe bouner Iranslorm. . ,N6 1 is equivalent lo two 

jf ' e k »* £ p(« - r) dr d« + jf ' ?(() e k " d; = <r{1) e 2 '"' df, (87) 
e k, » £ ^ <p(u - r) r,(T) dr du = ^ <r (0 e 21 "' dr. (88) 

Now, we observe dial (XX) does mil contain e and use it obtain an expression for i; as a 
function of (/,' and o. After that, we will view (87; as an expression for $ via <p, a and «/. 
l-'rom (XX) and the uniqueness of the f( airier transf. -nil, vae obtain 

j\(u -r)^(T)dr = io-Q), (89) 
for all « € [1, 2|. Introducing the new variable i via Ilia f. i inula 

» = u - 1, (90) 
we convert (89) into 

P( B+ ,- t )„( t )dr = Icr(!li). (91) 
which is a Vol terra equation ol the hrsi kind v, nli respect lo ,;; differentiating (91) with respect 

-<p(l) n(v) + f ¥ .>+l-r)/,(r)dr = \ a >( V -±±\, (92) 



ie is? i together with the uniqueness oi 



for all u e [0, 1]. 

I lie tollov.iiig theorem is a con sequence of the preceding one. 

Theorem 6.2. Suppose that a,tp ■ [-1.1] — C are a pair of e-'-fum-tio 
<p(-l) ^ 0, <p(l) # 0. e h a positive rail number ami i lie functions f p a, 
the formulae 



0 = J dt. 



Then there exist two functions ij. c : I 1.1] » C. such that 
for all x e R. with the functions q, r : [-1. 1] R defined by the formulae 



= J ^(t)e c "dt. 



Proof. IX-tining the function-. ;V. /' . /; , by ihe formulae 

/ + 0c) = jf orCOe^'df, (99) 

/ (.»:) = J aiO^'dl, (100) 

p + (x) = j^ <p{t)j""dt, (101) 

p_W = j'y(l)e™<At, (102) 

we observe that, for all jc s R 1 , 

/to = /+« + /-(*), (103) 

p(.v) = p + (.r) + p_(.v). (104) 

Due In theorem h.l . lliere e\isl such // , £ + , J , that 

/ t W = Pt W? t W + r t «, 1105) 

/_(.v) = P-W?-W + r_(.v), (106) 



with the functions <j+. <?_. r + . ;•_ defined by the formulae 
<?+« = J ij+Oe^'d/, 

9 00 = / >J 

r 0) = / I (/)e fa 'd/. 
Nov., defining <i by the formula 

p(x)q(x) = (p :v: I ;>.(v)i • (g_0) + $+(*)) 

and we dclme rl.vj h\ the obvious formula 

rO) = r_O0 + <" + O0 " (P-W« t W + /' "!</ (.11. 
1 he product ;M-V )<j (.v ) is given by the formula 



'••"«-™ -17* 



i/..in/j,ifie'"" +r 'd/dr; 



6.2. (Jtimirutwr unties firnit Ihe division ihe, Tern 

In much Ihe same way thai Ihe divWion theorem fur polynomials can he used to provide a 
constructive proof of ( iaussian quadratures, llieoreiu fi.2 provides a method ol constructing 
general i/ed (Iaussian quadratures for handd united li met ions. 'I lie method is us follows. 

'Jo construct a quadrature lor Junctions d a handv.adlli la . I'NYVI s corresponding to a 
baniKvidlh e are used, id bus the eigenvalues j.s, | a ndcigen Junctions j ,;. , | are in I his seel ion. as 
elsewhere in the paper, those corresponding to bandwidth ,-.) d he following theorem provides 
a bound of the error of a quadrature win >sc n ides are the n > its of [he mh prolate function f,„ 
when applied to a function ;" which satisfies the conditions of ihe division theorem, in terms 
of the norms of the quotient and remainder of / divided by fa : 

Theorem 6.3. Suppose thai jq.aj, r„ e R are ihe roots of fa. Let the numbers 

Y i w k fa(x k ) = j^faW&x, (115) 

for all j = 0,1 n - 1. Then for any function f : [-1,1] -»• C which satisfies the 

i-nnilitions of theorem 6.2, 

|£>*/0*)- fjWte^M-WriW + ll41-f>;MIV'jL'(2 + X><'*l). 

(116) 

where ihefttntiimis n. i" : I 1. I I <• € •«••• Jelhwii hi iluwm \2. 



Proot. Since f s iti fie I ulitioi theoren i rt tist functions q. r : [-1. 1] —>■ E 
defined by (97) and (98) such that 

llien, defining the error of integration E , for the function / by 



(118) 



E f = ^/Wcuj 
we have 

£/ = | E (V'„fe)9(^) + r(,: t )) - ^(V^W^W + rW)^! 

< | £tifc ifc (.*(-)« to) - j + l^u^fe) - ^ r(jc)<br|. 

(119) 

Since the nodes {.v t ) are the roots of t//„ , 

X! w k t» 0*t) * 0*) = 0- (120) 

B/ < ly'^W^Wdij + l^u^te)- ^ r (x)d*|. (121) 

y ^ v» (Jt )<? « = j t» m J i me™> dt dx 
= /i ?,(0 /i ^" Wei ™' <brd< 

= j W)K fn(t. 



h(t)dt. (122) 
.]ualu\ mid the lac Ilia the function lui> unit norm. \\c gel 



£ Ufcrfo) - J ^ r(x)dx = «*( J | ?(f )e u «' df) - j' ( £ ' f (Oe 1 ™' dx 

= /' eco( £ w * ci ""' - /' c '"' dA ) dt - (,24) 

Si ihsl liming (73 ) inlo i 1 24 i. and using the Caucln Scliw arl/ ncqi i . 1 1 il \ . v,o gel 
> 'I'l.ri v,; i 

< iifii-Ei^i-ii^i^-^+Ei^i)- (125) 



Combining '121). > 123 ! and ! 125). we get 

Ef « \K\-\\n\\ + \m-f^^j\-\\fj\t-(2 + ^\w k \j. (126) 

□ 

Remark 6.1. The use of theorem 6.6 lor the o instruction >f quadrature rules for band-limited 
fimctionv depends .in the fact thai the in this of flic handdiinited functions 4 and r in (117) are 
not large compared to the norm of f (both sets of norms being on [-%.. x-]>. Such estimates 
have been obtained for all n > 2c /x + 10 logic). The proofs are quite involved, and will 
be reported at a later dale. In this paper, we demonstrate the performance of the obtained 
quadrature kumulac nurncricalh i see section 8 below i. 

Remark 6.2. It is natural t< view 1 1 7 1 as an analogue for band-limited functions of the 
liuclid division theorem I or ]>ol ynoimals. I |.>\\c\cr. there are certain dil'lerenees. In particular, 
theorem 6.1 adrni Is extensions |,. Kind-limited fundi lis. f several variables, while the classical 
luiclid algorithm does not. Stich ex tensions ( la igel.lter v.utli se\eral applications i will be reported 



7. Interpolation via prolate spheroidal wavefunctions 



where <h u <t>2 4>» ■ [«. '•>] ^ C arc a fixed scqucn 

solve an i: / n linear system t.» determine the cocllicicn 
at the n inlerp(J;ilioii nodes, then use il27;U> evaluate f 



ii prolate spheroidal functions i// r) . 06 6',, i Tor a band limit < can provide a good 

regime where the accuracy is numerically useful, the error is of the same order of magnitude 

used as the interpolation fund ions in this procedure. I hey can be expected to yield an accurate 
interpolation scheme for baiiddnnilcd luiiclious. provided that the matrix to he inverted is 
well conditioned. The following theorem show s that, if the interpolation nodes are chosen to 
he quadrature nodes accurate up to tw ice the bandw dlh of interpolation, with the quadrature 
formula being accurate to more than twice as many digits as the interpolation formula is to be 
accurate to, llicn the matrix inverted m the procedure is close to being a scaled version of an 
orthogonal matrix. 

Theorem 7.1. Suppose the numbers u>,. "6 w» 6 R and x\,x-i x„ 6 Rare such 

I f L S 2ic "'dx-J2^ 2ic '" i \<^ d28) 



for all a e [- 1, 1], and for some c > 0. Let the matrix A be given by the formula 
VIM*,,) fix,) ■■■ f„-dx» 



Proof. Clearly 

where 6,, h the Kr^neeker del:..; lunelion. I Ivjnp (li> 

'* =a > - P' • (, ; /> ' ■ • 

= fy* - =— ^ J J (?) 1 (r) J2 u '<e~ k> " e io " dt dr. (133) 

Using (128), this becomes 

°» = a *-x^b£ £ f ^ (r) ■ ( £ ° c '" r d " ~ M+ r) ) d/ dr - 

(134) 

where / : |-2, 2] -» € is a funclii.n which satisfies ihc relation 

|/»(jc)| < f, (135) 
for all x <= [-2, 2]. Thus 

* = * - x-^ £ £ *- iw £ e_i "' e " dj d,dr 

+ - - J J ^j- ] (t)^ t - ] ij)f,(l + T)dldT. (136) 

Using (18), this becomes 

ejt=S ]k - j\j jW^ + - ^ ^ £fo i(r) £^ +r)d( dr. 

(137) 

Due In (he nrlhoniinnalily nf I lie fund ions | o .. |. llns becomes 



1 King I lie Cauchy-Schwartz inequality, this becomes 



; '' 4 L 1^-1 1 "^-'V /! I/ 1 , ^-iw^(' +T ) df | d 
H^I/££-~l^l 



needed fur a bandwidth of r and an 
for a bandwidth 2; and an accural-} 
for mlcrpuhilion of accuracy a . In ,.i 



le integrals of the interpolated 



8. Numerical results 

'Hie algorithms of sections 5 7 have been iniplcrncntcd in double precision i'4-hil boating 
point) arilliinefie, Willi the results shown in tables 1-4. lablcs I and 2 show the performance 
of quadrature nodes produceil by the schemes of sections 5 and 6. when used as quadrature 
nodes; 'tables 3 and 4 show tlieir performance when used as interpolation nodes. These are 
not actually the same sets of nodes; even with the bandwidth , ha- imcrp- 'Union being half of 
the bandwidth lor quadrature (as it is m the tables ). more nodes are needed to achieve a given 
accuracy of interpolation than are needed to achieve a given accuracy of quadrature, as can be 
seen by comparing the number of nodes (printed m the ciduinn labeled n in each table). The 
error figures in the tables are approximations of the maximum error of interpolation or of the 
quadrature, when applied to functions of the form cost,; i i anil sinio a i, with 0 < a < r; they 
were computed by measuring the error at a large number -1 points in a 1 for interpolation, in 
both n and or), including the extremes. The column labeled "Roots' o .mains the errors for the 
nodes produced by the scheme ol section fv the column labeled -Refined* contains the errors 
after those nodes, used as a starling point, have hecn run through the scheme of section 5. The 
variable ; which appears la the lablcs Is the requested accuracy, used to determine the number 
of nodes in the ways described in sections 5 and 7. 

Also tabulated arc the numbers of I egendre nodes required to achieve the same accuracy 
usi i i hi i n I n i i i I i u hem Si i c ( I J I c i lc regenerall 

known to he superior for interpolation, for thai case the numbers off liehyshev nodes required 
to achieve the same accuracy are also tabulated. 





l-lgurc 2 contains tlie iiJiix i i in i it i norm ofllic dumtiliw "I cadi prolate function i/.'d.rj. for 
c 200 and .re 1,1 1, as a function of j; also plotted, for comparison, is the maximum 





norm of the derivative of cticli normalized I .egendre polynomial /•,. i ) over the same range; 
and plotted below, on the same In iri/.. anal wale, are the absolute values of the eigenvalues kj. 
I lie graph shows that, lor this value of r. oompuiing the derivatives of a function given by a 
prolate series is a better-conditioned operation than computing the derivatives of a function 
given by a I .egendre series of the same ntiinher of terms, i Obviously, if the number of terms 
ean also be reduced, as in die situations oil tables I —I. there is a further improvement in ihe 
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condition number.! 'I he same general pattern of behaviour is exhibited for other values of c; 
sis r approaches zero (and I he prolate fund ions iij>j>riiiicl) the I .egendre polynomials), the value 
of/ at which ihe maximum norm of the derivative rises sharply also approaches zero (as is to 
he expected, since lor , i . the prolate lunclions reduce to I .egendre polj noinials). I'mally. 
tables 5 and 6 contain samples H'i|UudMlurv weights and nodes. 

Remark 8.1. In litis paper, detailed discussion of issues encountered in the implementation 
of numerical algorithms has been deliberately avoided, as well as any discussion of CPU 
lime requirements, mcmon requirements, etc. Thus, we limit ourselves to observing that all 
algorithms have been implemented in h'OR I RAN, that will) the exception of the procedure 
for the evaluation of PSWFs described in section 4 we have not designed or implemented any 
new or original numerical algorithms, and dial Ihe procedure of section 4 consists of applying 
standard tools of numerical analysts ..diagonali/ation of a tridiagonal matrix) to the well- 
known recursion (61). The resulting algorithm for the evaluation of PSWFs has the CPU time 
requirements proportional to r 2 , with a fairly large prop. aar main.) constant. Ihe procedure 



of|2|, wlien iipplia) li.llic system of l'uiicti.iii\ '!■ . '/'. ihvt. require order w' operations, 

also with a fairly large proportionality constant. On the other hand, the cost of finding all roots 
n of the function >,(•,. lyiny un lhc interval | - 1 . 1 1 is proportiiinal to /i. and the proportionality 
constant is not large. I lie largest i we have dealt will) in our experiments was about 6'dC'i 
with resulting quadratures having ahoiu Iv'KI nodes. In this regime, the construction of the 
quadrature I Mb nodes and weights) took several minutes on a 300 MHz Sun workstation; 
while Ihere are fairly ohvious ways lo reduce the cost of the calculation ihoth in terms of 
asymptotic ("PI ' lime requirements and m terms ol asvviukil proportionality constants i we 
have made no effort to do so. 

es presented in this section, ami 



quadratures obtained from the nodes of appropriately chosen profile lunclions: however, 
the tlijj'arw c between lhc numbers of nodes required by lhc two, approaches to obtain a 
prescribed precision is not large. When the nudes obtained via the two approaches arc 
used for the interpolation (as opposed lo lhc integration i of banddilililed functions, the 
performances of the two tire virtually identical. 

(2) For large c, the number of nodes required by a quadrature rale for the integration of 
banddimiied hi net ions w ill lhc band limit . is close lo -: the dependence on the required 
precision of integration is weak (as one would expect from theorem 2.? and subsequent 
developments). 

(3) The numbers of nodes required by our quadrature rales to integrate band-limited functions 

required by our inicrpohuioii formulae in order lo interpolate banddmiued functions is 
roughly nil limes less than the number of ( liebysliev (or Gaussian ) nodes. Again, the 
dependence of the required number . >f n. isles . at the accuracy requirements is weak. 

(4) The norm of Ihedillercnlialioii operator based on ournodes is of the order c 3 , as compared 
to the norm of the spectral differentiation operators obtained from classical polynomia 
expansions: ilns migln be useful in the design of spcclral on- pscudospcctrul i lech n iq lies. 




9. Miscellaneous properties 

PSWFs possess a rich sol of properties. vajtiiol\ rvscmhlinj: ihc properties of Bessel fnneiions. 
This section lists some of those properties. Some of the identities helow can he found 
in [5, 18,21]; others are easily derivable from the former. 
The identity 



(see section 5.1 litis a number of consequences which, while fairly obvious, seem worth 
recording, since similar properties of other special functions ha\c often heen found uselul. 



Differentiating < 1+ ) ) m times with respect to ,v and n times with respect to / yields the formula 

x-?«c k " = Q-j J^kj ,//"",,) o'/" ; (n. (141) 

for all .v. r e 1-1. 1|. Multiplying ! 140) by t~ u " J ' . and integrating with respect to r. com oris 

Sm( ' : ^~ U)) \ )_V to ^y(u). (142) 

faking die squared norm ufi 1411 1 ami imcgrating v.. ith respeel to a and r, yields the formula 

|]|A/ = 4; (143) 
combining this wilh (21) yields 

g,,, = |. (144) 

eir = E^^m- ( 145 » 

The identity 

= J e r " <lr s (l)6l (146) 

(see section 2.5) also lias a number of simple bin potentially useful consequences. 
Differentiating il k lames with respeel lo aa we gel 

Xjffix) = (ic) k j l e" , t k f j (f)it. (147) 
We next consider the integral 

f(x) = f(a,x) = f' ^-,/, j{ t)dt. (148) 
i hflcrcnlialina . I4X • wilh respeel lo i . v\e have 

^-f(a, x ) = icj^^-fj(l)6l. (149) 
Mulliplying i I4S) by ioo. and suhlracling il from ( 1 4'J i. we obtain 

±-f(a, x) - \caf(a, x) = ic£ c'"< it = \ckjfj(x). (150) 
In other words. /' satisfies the differenlial equaiion 

f{x) - icafXx) = ic>.j fj(x). (151) 
The standard 'variation of parameter' calculation provides the solution to (151 >: 

f(x) = ic\jj e- im( '-"Vo(Odi + /(0)e i '"™. (152) 



V = y o — (153) 
(i.e. V is the product of multiplication by 1 jie and differentiation ), we rewrite (147 > as 

V k (i/j)(x) = 7~ J ''t k " f/(Odt; (154) 
for an arbitrary polynomial i with real or complex coeliicients ). 

P(V)(tj)(x) - J' P(t)e cM f;(t)it. (155) 
By the same token, the function ip defined by the formula 

satisfies the differential equation 

P(D)(«W = ^,fcW. (157) 



Lemma 9.1. for any positive real i, hiteeer m > U o;«/.T € (-oo, +oo), 
O-* 2 ) +2) to - 2 (t + 1 )* +1) « + (/,„ - t (i + 1 ) - c 2 x 2 ) (x) 

-2c 2 kx tf _1) to - c 2 £(£ - 1) <tf~ 2! to = 0 (158) 
/or a» it ^ 2. Furthermore, 

(1 - ,v 2 ) VCM - 4a- <W + ( X „ - 2 - c 2 * 2 ) «// (.r) - 2c 2 .r ^to = 0. (159) 
In particular, 

-2 (* + 1) V'if +li (l ) + (Xm - k (k + 1) - c 2 ) ^'(1) 

-2 c 2 fc "(1) - c 2 k(k — I) ^* 2 ' ( 1) = 0 ( 1 60 ) 

for all k 2 2, and 

-2V/:(l) + (/,„-<t 2 )lMl) = 0, (161) 
-4^(1) + (x™ - 2 - c 2 ) <// (1) - 2c 2 ^(1) = 0. (162) 
Furthermore, for nil mleeer m -Delhi I -2. 

t/'if +2) (0) + (Xm -* (*+ 1)) >(0) - c 2 /fc(/fc - 1) ^" 2! (0) = 0. (163) 

For all odd rn, 

VC(0) + (z m :- 2)^(0) =0. (164) 
and for all even rn, 

V^(0) + X„V'»t(0) = 0. (165) 
Finally, for all integer m > 0 and k > 0, 

^(D^O, (166) 

tf£2i(°) = °. ( i67 > 

V'2 2 f +1) (0) = 0. (168) 



lMD = 0 
: integer m ^ 0 and observe that th 



for alii = 0, 1, 2 Due to the analyticity of •!:,, i.vi in the complex plane, this would imply 

V/,„W = 0 (171) 
foralLv e M 1 . □ 

The following is an immediate consequence of the identity (161) of lemma 9.1. 
Corollary 9.2. For all integer m, n > 0, 

ten • imd - *;,(D ■ iMD = (x. - xm) • ifr„(i) • *.d), (172) 

u7/c;v /,„, e t i «; _V 

'Ihcorcm 3.1, in seem in 3.1 . gives formulae for the entries oi matrices for differentiation 
of prolate series and lie multiplication of prolate series by ... Matrices for any combination 
of differentiation and multiplication by a polynomial can obviously be constructed from these 
two matrices; for instance, calling the differentiation matrix D, and the multiplicalion-by-.r 
matrix X, the matrix fur Inking Ihc second derivative ufa pod ate series, then multiplying it by 
5 -,v 2 , is equal to (5/ - X 2 )D 2 . 

In many cases, however, I lie re are simpler formulae for the entries of such mat rices, that is, 

lor inner products of i/ouv ] will] Us derivatives and Willi p. d\ la Is. I lie follow my theorems 

establish several such fonnulae. as well as a lev. ti>niualae f.-r inner products which do not 
involve ill o'.v ) 11 self bill only its denva lives. I heir proofs proceed along much ihc same lines 

as those of theorem 3.1 and are omitted. 

Theorem 9.3. ,V//o/>o.vc llnit • is real ami pasilivc, ami ilia; ihc uilcvcrs in ami n arc in -it- 



//m^«(mod 2), then 




j ^xf^(x)f m (x)dx = 0, 


(173) 


1 x 2 tZ<x)t m <x)ix =0, 


(174) 


J X 2 f' ll (.x)f' m (x)dx = () z 


(175) 




(176) 


f l x 2 f n (x)f m {x)dx = 0. 


(177) 



- (2 ^(1)^,(1) -«„„). 



// /;/ = m mod 2 i ami m * n, then 

j ^ x 2 iC u) i!,„ u) d x = 7-^77- « a ) ^(D-^CD V^ti)) 

= ^ K k (Xn - /,„) V'«d) <!'„,(i) - 3^r- V'*(i) tMD, 

J x 2 f' m (a) (x ) djr = 2 t' m (1) ifr„(l) + ^ 2> '" k {ir' m (1) V„ (1) - ^ (1) lfr m ( 

= 2 v/;; ( i ) .//,„ ( i ) + T^y- (*;; ( i ) ( i ) - v/-; ( i ) V'„ ( 1 x> 

= V'„,(DV' n ( 



j ^M*) i<l{x)<-\x = ; , 2 \, 2 cvCOD iMD - «/'„(!) ^(i)) 
= 77^7 (x„-/»)fc(i)^(i), 

j\r 2 v„ to f m (a i dA = - -§ Trr^ ( ^ ( 1 1 </'»> u ) - ^ u > ^ ( ' ) ) 

H-'frere Xm> X» e R are as defined in theorem 2.6. 



//>i is odd 1 and m is even, then 



1 lie abo\e theorems do mil. use much ol* Ihe detailed slruelure el i.lie mlcgral operators ol" 
which Ihe fund in us { i/; , 1 arc cij>cnfiim.'li»n<:. Tin is many of l he in genera I i/c o;isil> lo I lie ease 
ol an o|ieralor /, : f-|0. 1 | > /,-|0. 1 1 delined via the formula 

W)(x) = j\(xt)fU)dt. 1191) 

lor some function A' : |(>. 1 | > C; Ihe following llieoreni is an example of Ihis. 



Theorem 9.5. Let ) t tt operator L defined by i 19 1 >, that 

J K(xt)^(t)dt = ht,(x), (192) 
J K(xl)ii 2 {l)Al =X 2 f 2 (x)- (193) 

M J„ xf'(x) -AiUKlv 

provided liiai m'illier k, nor die deiinminalur of ///•' liyhl-hami \itlr <>l ( 194) is zero. 

'Ihc I'i'lli.win}! iliL-nrcni cM;ibJWhes ihc relation bclwccn ihc m-nii 'h'ench function >/• , on 
1 . 1 h which ru l.lns paper i> taken In he one ■ an J il> norm t>u i 

Theorem 9.6. Suppose that • is veal and positive, ami linn (he inteper a is non-nepaiixe. Then 

pf;(x)dx = \ (195) 

where p„ is given by (21). 

Proof. 

II W to = /I /' ^IV^ dr) W to 
= — y ^(0d* = — • 



J'* /to dJr = ^ / i M*) fMto. 



(196) 

Theorem 9.8. ''./'/'' d'.e >',,, 9 end p 'ad ", ,>„ <v ■.>>,".', t >, v,'i /h i' 

' '-ft. to. if — 1 < a < 1, 

W (197) 
if t > 1 or jc < -1. 



n of the operator Q c defined in (19), a 



p c 1 "' fjt)dt = j- J* cr*"(I I ^ Sm(C x ( * u " )} <M»)d«) dt (199) 

= JT / v ' m (!t) / ""^ -V "'" ei "' d ') dM - (200) 

Since the innermost integral is the orthogonal projection operator onto the space of functions 
of Hand limit r on i "v.. applied to the function e 1 ' il follows thai 

=i£ fc(a) ({r I ;L < i Jt o < r ';<-i}) d " 

Ji//™ ( " )e, ™ dM - ir -K- 1 ' (202 ) 

[ 0, if a: > i or x < -1, 

from which (197) foffows immediately. □ 

Hie following live theorems cslaHisli iTiniilac lor die derivatives id prolate functions 
and their associated eigenvalues with respect lo <■. Ponds "1 the lirst theorem can he found 
in [6,25]. 

Theorem 9.9. For all positive real • ami non-neyattve inteeei 

3K_ i ifllD-l 



^ = ^*i(». (204) 
The following theorem is an immediate consequence of the preceding one. 
Theorem 9.10. h'or nil positive real e ami non-ne'jalive integer hi. n. 

(¥) -X L 7^< 1 )-^"»- (205) 

Or) = 'jr \ _ (206> 

Theorem 9.11. Suppose that e is real ami posit ire. mill I lie iuiei;ers m. n are non-negative. If 



Proof. Since the norm of ijr„ on [- 1 . 1] remains constant as c varies, \jr„ must be orthogonal 
on [-1. 1] to its own derivative with respect to c. which immediately yields (208). To 
establish (207), we start with the identity 

X„f„(x) = j c WM f„(t)&t. (209) 
Differentiating <20V> with respect to r. we get 

£ fc (,) + ,„ ^ _ f (u^ M) + c - !M2) d , (2,0) 
Vlultiplviri" both Miles "l'i2IO> bv tfv w> arid inu-crulins with respect to .v. we have 

1 ./>. W ^ a -^/>*»«^^/>. l .)^* «.> 

which, using ( 1 /8l, wc rewrite as 

(K - K„) V',„(0 ^^d» = ^ (2^(1)^(1) - S m „). (212) 

Assuming lh;il ;,v == ,,. and dividing bv ).,, - /.„.. wc then get (2071. □ 
Theorem 9.12. Suppose that < /.» /<•«/ ami positive, ami the latepee ai is non-negative. Then 

(213) 



Proof. Due to theorem 2.6. 



-cV)V„W = 0. (214) 
= x™ + « and i//„(.t) = Vr m (jr) + S(x), this 



l + !(i)) = 0. (215) 
of the second order or greater that is. products 
and subtracting (214). wc get 
(1 - .v 2 ) S"(x) - 2xS'(x) + (Xm - c 2 x 2 )S(x) + (s - 2chx 2 )xlt m (x) = 0. (216) 
Let the sell-adjoint differential operator L be defined by the formula 

L(f)(x) = (1 - x 2 )f(x) - 2xf'(x) + (x,„ - cV)/fcr). (217) 



-. In addition, since /. is sdl-adjoinl. 

f\ L (X) fm(x)Ax = /, L «OMd*. (219) 

14), L(f m )(x) = 0 for all x e [-1. 1], so the integral (219) is zero. Thus (218) 

^-=2cj\ 2 r/ l 2 (x). (220) 



10. Generalizations and conclusions 



In this paper, we desi ul band-li ti 1 s based on the properties of 

PSWFs, and the connect! ms ol the latter « 1th certain fundamental integral operators (see (17) 
and (19) in section 2.5;. The quadratures are a surprisingly close analogue tor band-limited 
functions of Gaussian quadratures for polynomials, in that they have positive weights, arc 
optimal in the appropriately dehned sense, and their nodes, when used for approximation (as 
opposed In integration!, result in extremely efficient interpolation formulae. Thus, sections 5- 
1 of llns paper can ho Mewed as reproducing for baiid-limiled functions much of the standard 
polynomial-based appiexirnatnai theory I fur which sec, e.g., |26| I. Generally, there is a striking 
analogy bclwccn llv bdiiddiiml.cd bmcliuiis and polynomials. 

Obviously, there arc certain differences between Ihc resulting apparatus and the standard 
numerical analysis, lo slarl with, where the classical Icclnikjiics are optimal for polynomials, 
the approach of this paper is optimal lor band-limited functions: whenever the functions 
to be dealt with are naturally represented by trigonometric expansions on finite intervals, 
our quadrature and interpolation fumiulae lend lu be more efficient than those based on the 
polynomials. When Ihc functions lu be dcah with arc naturally rcprcscnicd by polynomials, the 
classical approach is more ellicienl; ||.»»e\cr. many physical pliouoinena involve band-limited 
funclions, and very few involve P'dy numials. 

Qualitatively, the quadrature . and rulerpolaiiun i nodes '.blamed in [his paper behave like 
a compromise bclwccn Ihc Gaussian nudes and the cquispaced ones: near the middle of the 
interval, lliey arc very nearly cquispaced. and near the ends, they concentrate somewhat, hut 
much less than Ihc Gaussian og-('hebyshe\ i nodes do. for large , . the distance between nodes 
near the ends of the interval is of the order of |/r' with the total number of nodes close 
to [ /jr. In contrast, the distance between the Gaussian nudes near the ends of the interval 
is of the order of l/ir, with n the total number ol nudes. A closely related phenomenon 
is the reduced norm of the differentiation operator based on the prolate expansions: for 
an n-point differentiation formula, the norm is of the order «••-, as opposed to n 2 for 
polynomial-based spectral differentiation. Ihus, PSWl s are likely lo be a belter tool for 
the design of spectral and pscudu-spcclral techniques than the orthogonal pely numials and 
related functions. 

Much of the analytical apparatus we use was developed more than 30 years ago 

(see 1 1 8, IU.2 1 .22] i; the f lame I importance ol ihese results ill certain areas of electrical 

engineering and physics has alsu been undersiuod for a long time. However, there appears 
to have been no prior attempt made lu view baad-linnlcd liiacliuns as a source of imiveriinl 
algorithms. Generally, there is a fairly limited amount of information in the literamre about 
llic I 'SWT's, especially when compared lo the wealth ol l ads on many other special functions. 
Section 9 of this paper is an attempl lo remedy ihis silualion lo a small degree. 

The apparatus built in this paper is a stricily one-dimensional one. Obviously, one can 
construct discretizations ri rectangles, cubes, etc. by using direct products ,.| one-dimensional 
grids; Ihc resulting numerical algorithms arc satisfactory but mil optimal, b'urthermorc, 
representation of band-hmhed functions un regtuus in higher iJiiiicnsiuns is ol both theoretical 
and engineering interest. Obvious applications include seismic da la collection and processing, 
antenna theory, NMK imaging, and many others. When the region of interest is a disc, 
most of the necessary analytical apparatus can be found in \11]. At the present time, we 
have constructed ami implemented somewhat rudimentary versions of the relevant numerical 
algorithms; we are conducting numerical experiments with these, and will report the results at 
a later dale. A much more difficult set of quest ions is presented by Ihc structure of band-limited 
functions on more general regions. 
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